Paleogene India-Eurasia collision constrained by observed plate rotation

Domain: earth_science
Published: 2023-11-10
Editorial outcome: minor_revision

Abstract
The Cenozoic India-Eurasia collision has had profound impacts on shaping the Tibetan plateau, but its early history remains controversial due to uneven availability of constraints. Recent plate reconstructions reveal two prominent counterclockwise rotation (azimuthal change) rate peaks of the Indian plate at 52-44 and 33-20 Ma, respectively, which could bear key information about this collision history. Using fully dynamic three-dimensional numerical modeling, we show that the first rotation rate peak reflected the initial diachronous collision from the western-central to eastern Indian front, and the second peak reflected the full collision leading to strong coupling between India and Eurasia. Further comparison with observation suggests that the initial and complete India-Eurasia collision likely occurred at 55±5 and 40±5 Ma, respectively, an inference consistent with key geological observations. We suggest that this collision history is instructive for studying the tectonic history of the Tibetan plateau and its surrounding areas.

Introduction
When India collided with Eurasia and how this process varied along the strike of the Indian front have been debated for decades1,2. Among the available constraints, the abrupt deceleration of the Indian plate’s drifting rate at ca. 50 Ma is often considered a direct response to the initial collision3–6. However, other geological proxies such as stratigraphy, sedimentology, metamorphism, and paleomagnetism along the suture often provide inconsistent ages of the initial collision2. A popular alternative explanation is diachronous collision7–12, but this model also bears large uncertainties with the proposed initial collision occurring in the west7,8, the central9–11, or even the east of the Indian front12. Better resolving the Paleogene collision history is critical for understanding the tectonic evolution of the Tibetan plateau and the associated geological and climatic implications, such as the plateau uplift11 and the monsoon evolution13.

The asymmetric distribution of land mass and seafloors within the Indian plate along the convergent boundary made the plate’s Cenozoic kinematic history highly sensitive to forces associated with its collision with Eurasia (Fig. 1 and Supplementary Fig. 1). The enormous collision-induced resistive force could not only reduce the speed of plate convergence, but also generates substantial resistive torques that may cause additional rotation of the incoming plate. Specifically, different suturing scenarios could generate distinct resistive torques, which may result in different patterns and rates of the Indian plate rotation. Observationally, some early paleomagnetic studies noticed significant Cenozoic rotation motions of the Indian plate1,14, but a continuous and quantitative description of this rotation history is lacking, so is the underlying geodynamic mechanism.

In this study, we analyze the Cenozoic rotation motion of the Indian plate based on multiple recent reconstructions15–18. We find two peaks in the rate of counterclockwise rotation during the Paleogene, which may bear new constraints on the process of the India-Eurasia collision. To quantitatively explore the mechanisms behind this rotation motion, we build fully dynamic three-dimensional subduction-collision models and test end-member scenarios of diachronous collision, including west-, middle- and east-bulging configurations of Greater India (the subducted northern extension of the Indian subcontinent)19,20. By comparing the observed rotation history with that predicted from the geodynamic models, we propose that the initial collision between India and Eurasia occurred around 55±5 Ma in western-central Greater India, and then propagated eastward, with the complete collision occurring around 40±5 Ma. We demonstrate that this collision process is consistent with a series of geological evidence in the tectonic domain.

Results
Cenozoic rotational motion of the Indian plate
To illustrate the Cenozoic motion of the Indian plate, we choose three reference points within the Indian sub-continent, which are reconstructed backward in time following a recent global plate reconstruction15. The trajectory and velocity of the central point are shown in Fig. 1a and c, respectively, which behave similarly to those reconstructed in previous studies3–6. In addition, we obtain the azimuths of the great circles connecting the northern and southern points, the variation of which describes the plate’s azimuthal rotation (Fig. 1b). We find that the Indian plate has rotated ca. 15.8° counterclockwise since the early Cenozoic, also consistent with previous estimates14. The rate of rotation exhibits two prominent peaks during the Cenozoic (Methods and Fig. 1c). The first peak occurred between ca. 52 Ma and ca. 44 Ma, accompanied by rapid plate deceleration, while the second peak occurred between ca. 33 Ma and ca. 20 Ma with an average plate speed close to that of the present. The rotation rate reached ca. 0.006°/km and ca. 0.010°/km during these two peaks, respectively, significantly faster than the mean background rotation rate of ca. 0.003°/km during Cenozoic (Supplementary Fig. 2d).

Besides the plate reconstruction of ref. 15, we further consider three other global reconstructions16–18 in estimating the Cenozoic motion of the Indian plate. These reconstructions differ not only in their relative plate circuits, but also in the absolute motions due to different reference frames (Supplementary Table 1). Interestingly, the resulting rotation patterns are generally similar, with slight differences in magnitude (Supplementary Fig. 2a–c). The first peak in the first model18 (Supplementary Fig. 2a) is not as robust as those in the other models, likely due to the discrepancy in the adopted absolute reference frames based on either the moving hotspots frame (e.g. ref. 15) or a slab-fitting scheme assuming vertically sinking slabs in the lower mantle21, where the latter was challenged by recent geodynamic modeling22. If the slab-fitting reference frame is replaced by the moving hotspots frame (e.g. ref. 17), the first peak in the first model will also appear similar to that derived from ref. 17. The additional rotation rate peak at ca. 60 Ma in some reconstructions (Supplementary Fig. 2a, b) is not considered in this study, since it is not observed in all reconstructions. The average of the rotation rate curves based on all four reconstructions exhibits two clear peaks (Supplementary Fig. 2d), where both the timing and magnitudes of these peaks remain almost the same as those in ref. 15. Consequently, we deem the two-peak pattern to be a reliable feature of the Cenozoic rotation history of the Indian plate.

Mechanisms for the variable rotation rates
To investigate the mechanisms that control the Cenozoic rotation of the Indian plate, we construct three-dimensional numerical models using the finite element code CitcomS following our recent effort23. A free-subduction system is designed where the Indian plate subducts beneath the Eurasian plate with an initial oceanic slab dipping into the mantle. To test the collision process, three different geometries of Greater India are considered, with a western bulge, middle bulge, and eastern bulge along its northern margin, respectively (Supplementary Fig. 3). Recent regional plate reconstructions24,25 suggest that the abutting oceanic lithosphere on the eastern side of the Indian continent is much wider than that on the west, and that both oceanic portions are separated from neighboring plates via weak spreading ridges or transform faults. This implies a continuous subduction history along the eastern side of the Indian plate during the Cenozoic, consistent with the recorded continuous arc magmatism in Sumatra26. A straight trench geometry (rotated to be E-W oriented along with the Indian plate, see Methods) is adopted at the model’s initial condition, as suggested from the along-strike paleolatitudes of the Lhasa terrane27,28, the restoration of intra-Asian shortening and extrusion19, and the existence of the linear high-velocity seismic anomaly in the lower mantle29. The overriding plate is fixed relative to the western boundary but free to deform on the eastern side, mimicking the rigid Eurasian continental blocks on the west and a relatively weak boundary on the east due to the Pacific and the Neo-Tethys subduction30,31. More details of the model setup are discussed in the Methods (Supplementary Figs. 3 and 4).

We simulate the collision process and monitor the rotation of the Indian plate. In the west-bulging model, Greater India initially collides with Eurasia on the western side and then sutures diachronously with Eurasia from west to east (Figs. 2a–d and 3b). The buoyant upper crust of the downgoing Greater India is scraped off along the suture, and its dense eclogitic lower crust and mantle lithosphere subduct continuously (Fig. 2e–h). In this model, the Indian plate experienced two periods of faster rotation, with the first reaching a peak rate of 0.0067°/km and the second of 0.0095°/km (Fig. 3a). Remarkably, both the timing and magnitudes of the two rotation rate peaks are consistent with the observation (Fig. 1c).

The two modeled rotation rate peaks postdate the respective initial collision and complete collision by about 8 Myrs (Fig. 2a–d), reflecting their formation mechanism. The reasoning can be explained by a simpler rotational torque analysis. We take the Indian plate as the frame of reference and the geometric center of the Indian plate as its center of mass, which is located next to the southeastern margin of the Indian continent. Assuming the along-strike variation of slab pull (per unit length along strike) is negligible, the collisional resistance will be the main force that influences the rotation of the Indian plate. During the initial collision, when the crusts of both continents first come into contact, the resistive force (F), assumed to be proportional to the length of the collisional front, and the rotational torque (T), a product of F and the length of the lever arm (L) relative to the center of mass, are minimal (Fig. 2a). This results in no noticeable extra rotation. Subsequently, the rapidly increasing length of the collisional front boosts F, whose effect dominates that of the decreasing L in determining T (Fig. 2b). The corresponding growing torque accelerates rotation, forming the first rotation rate peak. As it proceeds to a complete collision along the entire Greater Indian front, the decreasing lever arm L outweighs the increasing resistive force F (Fig. 2c), temporally reducing the rotation rate. This could also be verified through a simplified quantitative torque calculation (“Methods” section and Supplementary Fig. 5), where the torque T increases and reaches a peak first, followed by a decrease in T before the collision is complete. However, after the complete collision, F increases significantly due to the growing amount of accreted crust (thus gravitational potential energy or GPE), upper plate deformation, and inter-plate coupling (Fig. 2f–h and Supplementary Fig. 4b–d), causing T to enhance again (Fig. 2d), which generates the second peak of rotation rate.

For the middle-bulging model, the India-Eurasia collision occurs initially in the central part and then propagates both to the west and the east (Fig. 3c and Supplementary Fig. 6a–d). This model also exhibits two rotation rate peaks, with the first peak at a rate of 0.0049°/km and the second peak at 0.0080°/km (Fig. 3a). The two peaks postdate the initial and complete collision by about 6 Myrs, respectively. The amplitudes of the two peaks are slightly lower than those of the corresponding peaks in the west-bulging model. For the first peak, this is mostly because the length of the lever arm is shorter than that of the west-bulging model (Supplementary Fig. 6b). This explanation is also supported by the simplified torque calculation, in which the maximum torque (corresponding to the first peak) gets smaller when the bulge moves towards the east (Supplementary Fig. 5b). For the second peak, the difference in amplitude is likely caused by the E-W asymmetry in the distribution of resistive force and possibly slab pull force as well (Supplementary Fig. 6d). For the west-bulging model, the longer collision duration further west means a greater GPE and amount of upper plate deformation, thus a larger resistance west of the center of mass, relative to that in the middle-bulging case. In addition, it is expected that the slab on the east side of the west-bulging model is longer than that in the middle-bulging model following a larger accumulated rotation in the former, and the resulting greater slab pull further enhances the difference in their torques. However, given the uncertainties in both the reconstructions and the simulations (Fig. 3a and Supplementary Fig. 2), we suggest that the middle-bulging scenario is not necessarily inferior to the west-bulging one in terms of fitting the rotation records.

For a better comparison of different scenarios, we also test an east-bulging scenario where Greater India collides with Eurasia on the east first and then sutures diachronously towards the west (Fig. 3d and Supplementary Fig. 7a–d). This model exhibits only one rotation rate peak at a rate of 0.0061°/km, which is significantly lower than that observed (Fig. 3a). The first peak disappears in this model because the position of the initial collision is on the eastern side of the center of mass (Supplementary Fig. 7a), where the resulting clockwise torque does not allow an extra counterclockwise rotation. This clockwise torque is minor in magnitude and does not last long as the suturing propagates westward. The direction of the resultant resistive force is almost in line with the center of mass as the collision develops to Stage Two (Supplementary Fig. 7b), thus no apparent resistive torque exists to create any prominent extra rotation or the observed first rotation rate peak. The simplified quantitative calculation shows a similar tendency of torque variation when the bulge is located on the eastern side (Supplementary Fig. 5d). This model does generate the second observed peak, but with the lowest rate among all three cases (Fig. 3a), due to the similar reasons proposed in the middle-bulging case, including the E-W asymmetry in the distribution of the resistive and slab pull forces. Apparently, this model fits the reconstruction the worst. Therefore, an east-bulging Greater India is the least plausible scenario.

Implications for plate rotation and Paleogene India-Eurasia collision
Our proposed mechanism for the Cenozoic Indian abnormal rotation due to prominent torques could be further verified from the Mesozoic plate kinematics. For example, the Cretaceous Indian motion during its breakup from Gondwana32 also exhibits three rotation rate peaks, which occurred at ca. 140-121 Ma, ca. 100-83 Ma, and ca. 73-68 Ma, respectively (Supplementary Fig. 8a). The earliest peak corresponded well with the ca. 136-126 Ma separation between India and Australia-Antarctica32, likely triggered by the early eruption of the Kerguelen plume at the northeastern Greater India that formed the Comei large igneous province33 (LIP) at ca. 132 Ma. The second peak was consistent with the ca. 94-84 Ma breakup between India and Madagascar32 and the emplacement of the ca. 90 Ma Morondava LIP at Madagascar and southwestern India34. Similarly, the last peak was coeval with the ca. 71-62 Ma separation of Laxmi-Seychelles from India along its southwestern margin35 and the ca. 65 Ma eruption of the Deccan LIP34. The ca. 10 Myrs gap between the beginning of the rotation rate peaks and the LIPs eruption could be ascribed to the early effect of plume push before the plume head breaks through the lithosphere5. The localized plume push applied at the northeastern or southwestern margin of the Indian continent (Supplementary Fig. 8b) exerted large enough torque to have accelerated the counterclockwise rotation of the plate18. Mechanically, these plume-generated torques are similar to those due to the diachronous continental collision in driving the Indian rotation. The remarkable correlation between each rotation rate peak and an identifiable geodynamic torque since the early Cretaceous (Fig. 1 and Supplementary Fig. 8) strongly support the diagnostic implication of the observed Cenozoic Indian rotation on the history of India-Eurasia collision.

Given the good match of our modeled rotation history with that observed (Fig. 3a) and the fact that both the initial and complete collision stages predate the two rotation rate peaks by 5–10 Myrs (Figs. 2 and 3a and Supplementary Figs. 6 and 7), we propose that the India-Eurasia collision commenced at ca. 55±5 Ma along the western-central Tibetan margin and that the complete collision occurred at 40±5 Ma when the collision front expanded to the eastern margin (Figs. 3a and 4). These age ranges are consistent with various recent geological findings along the suture. In southern Tibet, sediment provenance changes have been reported in both deep-water strata of the northern Tethys Himalaya36,37 (Saga in Fig. 4a) and shallow-water strata of the southern Tethys Himalaya38 (Tingri in Fig. 4a), giving a collision age of ca. 59 Ma and ca. 51 Ma, respectively. In the western Himalaya, the same approach applied in shallow-water strata39 constrained the age of collision to be 54 Ma. Ultrahigh-pressure eclogite has also been reported in the western Himalaya40 with a peak-metamorphic age of ca. 47 Ma, suggesting a collision age at ca. 51 Ma. Further west in Pakistan, it has been proposed that the initial collision occurred at ca. 56-55 Ma based on provenance shift and tectonic responses41. Ultrahigh-pressure eclogite in this area with the age of peak-metamorphism at ca. 46 Ma has been reported, where the inferred collision age is similar to that in the western Himalaya42. In the eastern Himalaya, however, the collision age has not been well constrained due to the lack of stratigraphic record. High-pressure granulite from this region underwent prograde metamorphism from ca. 40 Ma, implying a similar collision age43. Further south in Myanmar, the collision of India-Burma with Eurasia commenced between the late Eocene and the early Oligocene based on ages of the unconformity44 and provenance changes45. In summary, although the geologically implied onset of the India-Eurasia collision varies along the suture, the ages generally fall between ca. 55±5 Ma and ca. 40±5 Ma, closely aligning with our model predictions (Figs. 3a and 4).

Various estimated ages of collision along the strike leave room for discussing the necessity and pattern of diachronous collision. Some studies believe that there is no diachroneity based on the comparison of stratigraphic records distributed in the central-western southern Himalaya2,46,47, but the collision in the eastern Himalaya and Myanmar is not constrained. Others support the existence of diachroneity, but with different preferences on the detailed collision process. Among them, some argue for an initial collision in the west, based on the earlier cessation of marine sedimentation and the earlier peak metamorphism in the west than those further east7,8. However, the inferred collision ages at different regions may not be intrinsically comparable due to the distinct characteristics of utilized rock samples43. Some recent studies based on multiple geological indicators along the entire suture proposed that the initial collision occurred in the center9,10. This logic, however, suffers from the fact that different indicators from the same region also produce variable ages of collision2. Besides, a recent study favored the earliest initial collision in the east based on the deformation feature in southeastern Tibet12. However, the deformation along the western-central Tibet was neglected in their discussion, casting doubts on the true along-strike variation implied from the analysis.

Utilizing three-dimensional fully dynamic numerical modeling, we find that the rotation patterns of the west-bulging and middle-bulging Greater India scenarios fit the reconstruction well (Fig. 3a). Based on the observed plate kinematics, our geodynamic analysis, and the existing geological evidence, we propose that the India-Eurasia collision could have either (1) initiated in the west and then propagated to the east, or (2) initiated in the center, then expanded to the western and eastern side (Fig. 4a). Synchronous initial collision in the central-western Indian front could also be a plausible option (label (3) in Fig. 4a), in which the central and western inner-margin areas (Gamba-Tingri and Zanskar in Fig. 4a) are considered to collide almost coevally, thus the central and western distal-margin areas (Saga and the western area with undiscovered deep-water strata in Fig. 4a) are believed to collide with Eurasia simultaneously2. According to our calculation (Supplementary Fig. 5), the first rotation rate peak will appear after suturing around or more than half of the whole Indian front. The distance between Gamba-Tingri and Zanskar (thus between Saga and the undiscovered area) is around one-third of the whole front (Fig. 4a), making it possible to allow two rotation rate peaks. The time interval between the initial collision and the first peak in this scenario will be shorter than in the bulging ones, but still lies in the range of uncertainty and thus implies a similar age of initial collision. These scenarios with western-to-middle-bulging Greater India are compatible with the inferred collision ages discussed above and the widely accepted collisional process suggested by earlier studies2,7–11,46,47. Several recent paleomagnetic studies sampling along the Tethys Himalaya delineate the shape of Greater India and suggest that its western-central margin extended farther north than its eastern margin20,28, further supporting our geodynamically inferred diachronous collision.

While the first rotation rate peak correlates well with syn-collisional processes, the second rotation rate peak during ca. 33-20 Ma (Fig. 1c) might correspond to a special post-collisional stage. All three numerical models exhibit rotation rate peaks during this period (Fig. 3a), which is mainly ascribed to the increasing coupling between the two plates after the final suturing, consistent with the gradual uplift of Himalaya and growing GPE during Oligocene11. We notice that this strong plate coupling corresponds well with enhanced tectonic deformation in southeastern Tibet48,49. The two major continental blocks in this region, Indochina and Sibumasu, underwent extrusion of hundreds of kilometers following a clockwise rotation of dozens of degrees48 (Fig. 4b). The ages of slip along the shear zones between blocks have been constrained to range from early Oligocene to early Miocene48,49, coeval with the post-collisional rotation rate peak of this study. In addition, the seafloor spreading of the South China Sea50 initiated at ca. 33 Ma and terminated at ca. 16 Ma, which might also be related to the contemporaneous extrusion of these blocks30,31 (Fig. 4b). We suggest that the second peak of the Indian plate rotation rate and the intense deformation within southeastern Tibet during ca. 33-20 Ma both reflected strong plate coupling after the complete India-Eurasia collision, which profoundly influenced the evolution of the Tibetan plateau.

There is a possibility that the India-Eurasia collision was more complex than a single-collision scenario51. Some of the proposed collision hypotheses can be compared against the rotation history of the Indian plate. Kohistan-Ladakh and West Burma are regarded in several reconstructions52,53 as oceanic arcs, which could have formed in a double-subduction system situated within the Neo-Tethys Ocean during the Cretaceous24. This system is believed to exert a larger slab pull on the adjacent plate than a normal single-subduction system54,55, such as that located at Sumatra in the east. If true, there would have been a torque resulting in an additional clockwise rotation of the Indian plate until the arc-continent collision eliminated the asymmetry in slab pull. In this case, a continuously increasing rotation rate (assuming negative for clockwise rotation and positive for counterclockwise rotation) is expected from the weakening of this clockwise torque, followed by the strengthening in counterclockwise rotation induced by continent-continent collision. The Cenozoic double-peak rotation rate pattern (Fig. 1c) seems inconsistent with the scenario where arc-continent collision happened after ca. 60 Ma (e.g. refs. 56). However, arc-continent collision prior to ca. 60 Ma (e.g. refs. 57) is a plausible scenario, although its rotation rate effect might have been obscured by the influence of plume push in the late Cretaceous (Supplementary Fig. 8). Besides the arc-continent collision, several other configurations of Greater India have been proposed, which generally assumed that this region had undergone significant extension before the collision19,58. We suggest this would not change the major modeling results, since the stretched continental lithosphere or the new oceanic lithosphere is still more buoyant than the surrounding old oceanic lithosphere58, thus more resistive upon collision. Therefore, we speculate that Greater India with a normal or stretched lithosphere are both compatible with the observed rotation features.

To conclude, the Cenozoic rotation motion of the Indian plate could be a novel constraint on the otherwise elusive history of the Paleogene India-Eurasia collision. Among the two prominent counterclockwise rotation rate peaks, the one during ca. 52-44 Ma could reflect the diachronous collision that initiated in the western-central front at 55±5 Ma and then sutured towards the east until 40±5 Ma. The second rotation rate peak during ca. 33-20 Ma should have recorded the formation of strong plate coupling after the final suturing, which drove the deformation of southeastern Tibet. This model of Paleogene India-Eurasia collision is consistent with various geological observations across the tectonic domain. These results suggest that plate rotation history has strong diagnostic implications on abnormal plate kinematic events and the underlying forces and torques that drive these events.

Methods
Rotation rates calculations
We calculate the Indian plate rotation rates by measuring the azimuthal changes of the great circle that connects the northern and southern points normalized by the distances that the plate traveled (blue lines in Fig. 1c). This definition of rotation rate is warranted for multiple reasons. First, it can mitigate the effect of the fast time-varying velocity of the India plate (dotted green lines in Fig. 1c), which has an inherent azimuth variation component, except in special circumstances, such as when the Euler pole is located at the North/South pole, or at the equator 90° of longitude away from the central reference point59. Second, the Cenozoic Indian plate rotation concerns mostly the dynamic effects of the east-west asymmetry of the downgoing plate at different stages of subduction, which is irrelevant to how fast the plate traverses (Torque calculations in Methods).

Numerical model
We use the 3D spherical finite element code CitcomS60 to simulate plate subduction and continental collision. The code solves for thermal-chemical convection governed by the conservation of mass, momentum, and energy, with the assumption that the mantle is incompressible and satisfies the Boussinesq approximation.

The model covers a domain of 72°×82°×2890 km in longitude×latitude×depth, which is discretized into 768×768×128 elements. The size of elements varies spatially, with the finest resolution of 10×8×8 km occurring in the vicinity of the subduction and collision zones. The model’s inner domain, with a size of 50°×72°×2890 km, covers parts of the eastern Asian continent, the Indian continent, the Neo-Tethys Ocean, and the Indian Ocean before the continental collision. All features are extracted from a recent plate reconstruction15 with necessary simplifications for the initial condition (Supplementary Fig. 1), and are rotated counter-clockwise to ensure the model’s N-S orientation is parallel to the direction of plate convergence whose actual azimuth falls between ca. 20°−40° in plate reconstructions15–18. The azimuth of the convergence boundary was ca. 290°−310° estimated from the geological proxies19,27–29, largely perpendicular to the direction of convergence. Thus the trench is correspondingly rotated to E-W direction in the models. A weak region beyond the inner domain is applied to provide enough room for Indian rotation and to avoid artificial return flow from side walls. All boundaries of the model are free to slip which allows for natural subduction and realistic rotation (Supplementary Fig. 3). The basic model setup is similar to that of ref. 23.

Three-dimensional depth-, temperature-, strain-rate-, and composition-dependent viscosity are considered in the model. A four-layer background viscosity profile is assumed here, including the lithosphere, asthenosphere, transition zone, and lower mantle. Their absolute viscosity values in the models are on the order of 1022−1023, 1019−1020, 1019−1020, and 1022 Pa·s, respectively. Thin weak layers are applied at oceanic subduction interface to mimick dehydration weakening23, whose absolute viscosity is set to 8×1018 Pa·s (Supplementary Fig. 4). The generalized viscosity law used here is

η(T,r) = A(r)ηc ϵ_II^(1/n - 1) exp [Ea/(nR(T + Toff)) - Ea/(nR(Tm + Toff))] (1)

where η is viscosity, T is temperature, r is radius, A is the pre-factor of viscosity, ηc is the pre-factor of composition-dependence, ϵ_II is the second invariant of the strain rate tensor, n is the creep exponent, which equals 1.0 for diffusion creep and 3.5 for dislocation creep, Ea is the activation energy, Toff is the temperature offset, Tm is the mantle temperature and R is the gas constant. The viscosity of continental lithosphere is controlled by composition tracers considering the complexity of continental geotherm and composition, and its temperature effect is not considered for simplicity. All the depth-dependent rheology parameters are detailed in Supplementary Table 2. We use a relatively low activation energy for diffusion creep compared to values suggested by laboratory experiments61,62, to account for various slab weakening mechanisms such as slab faulting, hydration, and grain-size damage63. The final effective viscosity is controlled by the combination of diffusion creep, dislocation creep, and pseudo-plasticity

ηeff = min(σy/ϵ_II, ηdif + ηdis) (2)

where the diffusion creep, ηdif, is applied in the whole domain, while the dislocation creep, ηdis, is only considered in the domain shallower than 660 km. σy is the yield stress following the Drucker-Prager yield criterion64,65

σy = μP + C (3)

μ = μ0(1 - min(1, εp/εf)) (4)

C = C0 + (Cf - C0)(1 - min(1, εp/εf)) (5)

where μ is the coefficient of friction, μ0 is the initial coefficient of friction, C is the cohesion, C0 is the initial cohesion, Cf is the minimum cohesion, εf is the reference plastic strain, εp is the accumulated plastic strain after the yield stress is reached. The material gets weaker as plastic strain accumulates. The lateral composition-dependent rheology parameters are detailed in Supplementary Table 3.

The composition field is defined using Lagrangian tracers (Fig. 2, ref. 66). The continental lithosphere has a three-layer structure which consists of the upper crust, the lower crust, and the lithospheric mantle. The Eurasian continent is further divided into the weaker non-cratonic and stronger cratonic regions, with the latter including the Karakum, Tarim, Qaidam, and Sichuan basins. The oceanic lithosphere has a two-layer structure that consists of the crust and the lithospheric mantle. For convenience, the initial temperature fields of the continental and oceanic lithosphere are both defined using the half-space cooling model, assuming a mantle potential temperature of 1300°C. The initial age is 50 Myrs for the Indian oceanic plate, 100 Myrs for the Eurasian non-cratonic regions, 150 Myrs for the Indian continent, and 200 Myrs for the Eurasian cratonic regions. Temperature-, composition-, and phase-change-dependent density anomaly is considered in the models. The density anomaly used here is defined as:

δρ = -α0ρ0ΔT + δρphΓ + δρch (6)

where α0 is thermal expansivity, which equals 3.0×10−5 /K, ρ0 is reference density, which equals 3340 kg/m3, ΔT is temperature difference, δρph is the density jump across the phase change, Γ is the phase function67, δρch is the chemical density difference between the compositions. The detailed thickness and chemical density setups for all compositional layers are listed in Supplementary Table 4.

Torque calculations
We perform a straightforward calculation on the variation of the collision-related resistive torque (Supplementary Fig. 5). We take the Indian plate as the frame of reference. The width of the Indian continent is s, which equals 3500 km in the calculation. The center of mass of the Indian plate initially lies s/3 away from the eastern boundary (represented by parameter b in Supplementary Fig. 5a). Considering the counterclockwise rotation of the Indian plate (e.g. Figure 2) and asymmetric trench geometry due to the continental indentation on the west30,31 and trench retreat on the east68, the unsubducted portion of the plate on the west is longer than on the east. The center of mass moves westward linearly as the collision goes on, with the maximum offset reaching 500 km. The southward movement of the mass center of the remaining Indian plate is not considered in the calculations, since this does not affect the length of the lever arm and the magnitude of the resistive torque. The three parameters above are estimated according to the recent reconstruction15. The location of the bulge varies from 0 to s away from the western boundary (represented by parameter a in Supplementary Fig. 5a). The resistive torques are calculated via the contact width of the continents, thus the effects of the plate’s time-varying transverse motion and the latitude-varying locations of the bulge are eliminated. We take the resistive force along the collisional front69 as 8.0×109 N/km, and assume the force is proportional to the contact width.

When the bulge lies at the west side of the collision front (a < s/2), the torque increases quickly after the initial collision (contact width = 0 in Supplementary Fig. 5b) and reaches a maximum peak, then decreases gradually (blue, orange, green and dark red curves in Supplementary Fig. 5b). This scenario is similar to the first rotation peak in the west-bulging model (Fig. 3a, b). In the middle-bulging model (a = s/2, purple curve in Supplementary Fig. 5b), due to westward movement of the mass center and associated shortening of the lever arm, its behavior is similar to the west-bulging model but with a lower first peak (Fig. 3a, c). When the bulge lies at the east side of the collision front (a > s/2), the center of mass in the model moves eastward first and then moves westward. Thus the total displacement of the mass center is very minor (Supplementary Fig. 5c). The corresponding calculation is shown in Supplementary Fig. 5d. The torque decreases slightly and reaches a minimum, then it increases quickly (pink, gray, and light green curves in Supplementary Fig. 5d). This scenario is similar to the East-bulging model that does not exhibit the first rotation peak (Fig. 3a, d). We emphasize this torque calculation is a first-order simplification of model results. The torque of east-bulging scenario with a fixed mass center during the full collision (e.g. gray curve in Supplementary Fig. 5d) reaches a larger value than that of west-bulging scenario with a moving mass center (e.g. orange curve in Supplementary Fig. 5b), which seems to be contradictory with the model results (Stage Three in Fig. 3a). This could be ascribed to the increasing resistive force on the eastern portion of the east-bulging model as collision goes on, which would weaken the final counterclockwise rotation.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is an interesting and innovative paper which links periods of enhanced rotation of the India plate with the India-Asia collision history. Overall, the proposed scenarion makes sense, I think, and, to my knowledge, this link has not been proposed before. The modelling is state-of-the art, and I also like the simple torque computations which help understanding what is going on, and the paper is quite well-written.\nI think there are two main points (and a couple of minor ones) though, where the paper could be improved I think.\nMy first issue is that I think the result, and in particular the first rotation peak, is not very robust. It occurs in the Müller et al. models which are, I suppose, related as they are both with Müller als first author but cannot clearly be seen in the others (not at all in van Hinsbergen, only very vaguely in Torsvik). Because there are two out of four models from Müller, it can also be seen in the average, but I think this doesn't mean anything. What I am wondering is whether, apart from comparing different models, it is also possible to do a formal error analysis? I mean, plate motions and rotations usually come with uncertainty ellipses, so would it be possible to compute whether these rotation peaks are significant, considering uncertainties?\nMy other issue is that some of the statements and figures in the paper appear actually contradicting each other. Maybe I misunderstand something here, but it requires clarification: In Figure 3a, all three model curves show counterclockwise (positive) rotation all the time. But in the other three panels, there is only counterclockwise rotation seen in b (west-bulging model), whereas c (middle-bulging model) shows more or less no rotation at all, and d (east-bulging model) actually shows clockwise rotation. And the middle-bulging model presumably corresponds to the a=s/2 case in Extended Data Figure 5. But in contrast to Fig. 3a, this shows no peak but just a continuous increase. And also the East-bulging model looks different again in Extendend Data Figure 5: Here the torque is first negative, then becomes positive, whereas in 3a it is always positive and in 3c always clockwise (i.e. negative). So, I am confused.\nFurther comments:\nline 114: But, North-South convergence, as described in the Methods (l. 499/500) is not identical to E-W boundary, unless the convergence is orthogonal to the boundary. Is that the case?\nline 238: Isn't the pre-dating by 5-10 Myrs best seen in Fig. 3a?\nline 265: Why is Sibumasu Block abbreviated with \"SP\"? What does the \"P\" stand for?\nline 482: \"as long as the Euler Poles do not overlap with the North Pole\" - I am not sure what you mean with this. I mean, also if there is an Euler pole at the equator 90° away from the continent, this corresponds to pure northward motion with no rotation.\nFor the absolute viscosity values in line 508, do you mean in reality or in your model? I am also wondering whether you need to specify a higher viscosity in the lithosphere; doesn't that come in automatically with the higher temperature, so you don't need to explicitly prescribe a higher viscosity (because otherwise, the same effect is considered twice)?\nIn your viscosity law in line 512, everything seems to be dimensionless, so how about the second invariant of the strain rate tensor, this should normally be in units of 1/s; how do you get this dimensionless?\nline 571: I think the sentence \"The width ...\" can be improved. I prefer how it is said on line 109/110 \"oceanic lithosphere on the eastern side of the Indian continent is much wider than that on the west\"\nExtended Data Figure 5: Why is there a \"1e16\" on top of panel b? Why you don't include curves for a=0 and a=1?\nExtended Data Table 3: So here you don't distinguish between Upper Crust, Lower Crust and Lithospheric Mantle? An the distinction between Asian Non-Cratonic and Cratonic Lithoshere is for all three (Upper Crust, Lower Crust and Lithospheric Mantle)?\nExtended Data Table 4: Does \"Neutral\" mean it is the same as the underlying mantle, apart from Temperature dependence? And what is actually the temperature dependence of density (i.e. thermal expansivity)? I don't see it anywhere specified.\nMinor commenst:\nline 33: \"has been debated\"\nline 294: \"could have either\"", "reviewer_2": "This is an interesting paper bringing some new constraints to understanding the history and potential driving forces of the India-Eurasia collision. Using the connection between changes in the azimuth of the India-Eurasia collision to constrain geodynamic models is a novel way of looking at this problem. I enjoyed reading the paper, and suggest that it should be published with some relatively minor corrections. The English language needs a bit of improvement throughout the paper. I have made some comments along these lines.\nIn terms of the interpretation of the results, I have this overall comment: The authors are following the premise that the time-dependent changes in rotation of India relative to Eurasia can only be explained by a successive continent-continent collision which needs to start around 55Ma or earlier, necessitating an at least 2000 km wide thinned continental margin off northwestern India (the details of the inferred geometry are difficult to extract from the paper).\nThe first question one needs to ask oneself is: Is this geodynamically reasonable or possible? The widest continental margin in existence today is about 400 km in width. There is a simple reason for this. It is not rheologically possible to create a much wider passive margin, as the stretched lithosphere will simply break, producing seafloor spreading, preventing excessive stretching from occurring. Stretching factors larger than four lead to SFS (see classical papers by Dan McKenzie). Continental crust is too weak to be stretched beyond this limit. So how exactly do you propose to generate a 2000 km wide passive continental margin to start with, 5 times wider than any margin in existence? In my view this essentially falls in the land of geodynamic phantasy, so I'd like to see this justified in the paper (personally, I don't think this is possible, but I am looking forward to see the explanation you come up with!). The subduction of this much continental crust also presents other problems, as is discussed, for instance, in van Hinsbergen DJ, Lippert PC, and Huang W. Unfeasible subduction? Nature Geoscience. 2017 Dec;10:878-9.\nAlong the same lines, I am surprised that the recent paper by Pusok and Stegman (The convergence history of India-Eurasia records multiple subduction dynamics processes, Sci Adv 2020) is not cited at all. Essentially it provides some ideas for an alternative way to create a two-stage collision history, with the first collision reflecting the collision of an island arc with Eurasia followed later by a continent-continent collision. This scenario was also proposed in a number of other papers. This should at least be discussed as a potential alternative way to explain the kinematic observations.\nSome detailed comments:\nLine 16:\nThe Cenozoic India-Eurasia collision has profound impacts on the formation of the\nmodern Tibetan plateau.\nshould be\nThe Cenozoic India-Eurasia collision has had profound impacts on the formation of the\nmodern Tibetan plateau.\nLine 136, and many subsequent instances:\nI do not understand why in \" West-bulging model\", \"Middle-bulging model\" etc \"west\" and \"middle\" is\ncapitalized. It should be lower-case.\nLine 145:\nThe two modeled rotation peaks postdate the respective initial collision and complete\ncollision by about 8 Myrs (Fig. 2a-d), which bears their formation mechanism.\nshould be\nThe two modeled rotation peaks postdate the respective initial collision and complete\ncollision by about 8 Myrs (Fig. 2a-d), reflecting their formation mechanism.\nLine 150: neglectable\nshould be\nnegligible\nLine 154: This results in no noticeable rotation.\nThe use of the word rotation is ambiguous here, as well as through much of the rest of the paper. The\nissue is that all plates rotate about Euler poles at all times, ie plates are always in rotation. What you\nmean here is rotating about a pole located proximal to the Indian continent. You need to think about\nhow to do that better, ie you cannot write that at some point India was not rotation, and then it started\nto rotate, seeing that it is ALWAYS rotating.\nLine 166: Similarly, a \"peak in rotation\" is not the correct term here. You actually mean a peak in\nrotation rate. The same issue occurs in other places in the text.\nLine 261: Some of the letter labels overlying \"Greater India\" are barely visible, especially those\nfrequently referred to in the text like the \"Saga\" and \"Tingri\" areas."}, "author_response": "Response to Reviewers’ Comments\nIn this response, we provide point-by-point responses (in blue) to all the comments (in black).\nThe revised manuscript is marked in red and green displaying the revisions.\nReviewer #1 (Remarks to the Author):\nThis is an interesting and innovative paper which links periods of enhanced rotation of the\nIndia plate with the India-Asia collision history. Overall, the proposed scenarion makes sense,\nI think, and, to my knowledge, this link has not been proposed before. The modelling is state-\nof-the art, and I also like the simple torque computations which help understanding what is\ngoing on, and the paper is quite well-written.\nI think there are two main points (and a couple of minor ones) though, where the paper could\nbe improved I think.\nThank you very much for taking the time reviewing our manuscript. We appreciate your\npositive comments and valuable suggestions on our paper that have greatly improved our\nmanuscript.\n(1) My first issue is that I think the result, and in particular the first rotation peak, is not very\nrobust. It occurs in the Müller et al. models which are, I suppose, related as they are both with\nMüller as first author but cannot clearly be seen in the others (not at all in van Hinsbergen,\nonly very vaguely in Torsvik). Because there are two out of four models from Müller, it can\nalso be seen in the average, but I think this doesn't mean anything.\nThis is a good point. We have re-examined the rotation parameters of the four plate\nreconstructions and tested the effects of different relative plate circuits and absolute plate motions\n(Extended Data Table 1). As for the relative plate circuits used by van Hinsbergen et al. (2021)\nand Torsvik et al. (2019), they are actually based on the same reference papers for relative plate\nmotions, with few differences in geologic time scale. The reason why their rotation features\nbehave differently is mainly due to the different absolute motion frames they use. For example,\nvan Hinsbergen et al. (2021) used the slab-fitting absolute frame that reconstructs plate motions\nassuming slab remnants sink vertically in the lower mantle (van der Meer et al., 2010). This\nassumption is highly controversial considering the enormous lateral slab displacement caused by\nhorizontal mantle flow (Peng & Liu, 2022). If we combine the plate circuits from van Hinsbergen\net al. (2021) with the absolute motion frame used by Torsvik et al. (2019), the \"hybrid\"\nreconstruction is similar to those derived from Torsvik et al. (2019) (Fig. R1), thus verifying that\nthe two relative motions in van Hinsbergen et al. (2021) and Torsvik et al. (2019) are indeed\nsimilar, both showing a peak during 52-44 Ma, although not as significant as in Müller et al.\n(2016) and Müller et al. (2019). We prefer the reconstructions using absolute frames based on\nmoving hot spots and various other constraints, such as those proposed by Torsvik et al. (2008),\nDoubrovine et al. (2012), and Tetley et al. (2019) that are used by the other three plate\nreconstructions (Müller et al., 2016; 2019; Torsvik et al., 2019). Here, we want to emphasize that\nexcept the controversial slab-fitting absolute frame, combining\nother absolute frames (e.g. the three mentioned above) and different relative plate motions\ndoes generate the peak during 52-44 Ma. We show van Hinsbergen et al. (2021) for\ncompleteness. Consequently, the average of the four reconstructions with different relative\nplate circuits and absolute frames effectively weighs down the contribution from the slab-\nfitting absolute frame, and thus should better represent the preferred plate motion. We have\nadded some extra descriptions in the main text (lines 94-98).\nFig. R1. The rotation rate and velocity derived from the \"hybrid\" reconstruction (van\nHinsbergen et al. (2021) relative motion + Torsvik et al. (2019) absolute motion) (left) and the\nTorsvik et al. (2019) reconstruction (right).\nWhat I am wondering is whether, apart from comparing different models, it is also possible to\ndo a formal error analysis? I mean, plate motions and rotations usually come with uncertainty\nellipses, so would it be possible to compute whether these rotation peaks are significant,\nconsidering uncertainties?\nWe appreciate your constructive suggestion. We have collected all available covariance\nmatrices of the finite rotations used in the four global reconstructions, which denote\nuncertainties in reconstructions (Kirkwood et al., 1999). But the covariance matrices of absolute\nmotions are lacking (e.g. Tetley et al. (2019) for Müller et al. (2019), Torsvik et al. (2008) for\nMüller et al. (2016), van der Meer et al. (2010) for van Hinsbergen et al. (2021)), which\nprohibits the evaluation of absolute motions of the Indian plate. Therefore, we are unable to\nconvert the segmented rotation uncertainties to a composite uncertainty or to get their\nconfidence region. The calculation for the uncertainties of rotation rates have not been studied\nbefore, which is beyond our capabilities and scope of this study. Besides, as suggested by White\n& Lister (2012), the true uncertainty of reconstruction involves many more factors than are\ncurrently taken into account, such as the deformation within the plates and the uncertainties\nwhen distinguishing magnetic anomalies. We will relay this effort to future research.\n(2) My other issue is that some of the statements and figures in the paper appear actually\ncontradicting each other. Maybe I misunderstand something here, but it requires clarification: In\nFigure 3a, all three model curves show counterclockwise (positive) rotation all the time. But\nin the other three panels, there is only counterclockwise rotation seen in b (west-bulging\nmodel), whereas c (middle-bulging model) shows more or less no rotation at all, and d (east-\nbulging model) actually shows clockwise rotation.\nThanks for pointing this out. We agree that this figure is a bit misleading. We have redrawn\nFigure 3 to make it clearer. We have added the subducted portion of Greater India at the\nstages two to four with dashed lines, which better exhibits their trajectory after collision. The\nnortheastern edge of the plate moves faster than the northwestern edge in all three models,\nwhich denotes their continuous counterclockwise rotation.\nAnd the middle-bulging model presumably corresponds to the a=s/2 case in Extended Data\nFigure 5. But in contrast to Fig. 3a, this shows no peak but just a continuous increase. And\nalso the East-bulging model looks different again in Extendend Data Figure 5: Here the torque\nis first negative, then becomes positive, whereas in 3a it is always positive and in 3c always\nclockwise (i.e. negative). So, I am confused.\nThis is a very good point. The simple torque analysis is really a simplification of the complex\ncollision process simulated by our models. We have updated the torque calculation (Extended\nData Fig. 5), which should be more realistic and reasonable than the previous analysis (lines\n164-166, 179-182, 204-205, 610-614, 626-643). We notice that the mass center of the Indian\nplate would move westward after the collision both in the numerical models and in the plate\nreconstructions. Because of the extra counterclockwise rotation since the collision, the Indian\nplate remaining on the surface would become wider on the western side than on the eastern\nside in the numerical models (Fig. 2). Thus the center of mass should move westward along\nwith the collision process. This tendency is also reflected in plate reconstructions (Gibbons et\nal., 2015; Zahirovic et al., 2016). The eastern boundary of the Indian plate, the Wharton Ridge,\nintersected with Eurasia from Java and eastern Sumatra to western Sumatra during 55-40 Ma.\nThe displacement of the ridge is estimated to be ca. 1000 km, which corresponded to the ca.\n500 km westward movement of the center of mass. Besides, the trench advance in the west\ncaused by continent indentation (Tapponnier et al., 1982; Jiao et al., 2023) and trench retreat in\nthe east due to back-arc extension (Li et al., 2013; Sternai et al., 2014; Capitanio et al., 2015)\nwould also facilitate the asymmetry of the remaining plate [Redacted].\n[Redacted]\nThen we add this feature into the torque calculation, in which the center of mass moves\nwestward continuously during the collision process (Extended Data Fig. 5, Fig. R3a,b). The\nresults are more consistent with the modeling results, with both the west-bulging scenarios\n(a<s/2) and the middle-bulging scenario (a=s/2) exhibiting positive peaks of resistive torques,\nwhich correspond to the first peak of the rotation rate (Fig. 3a). Most east-bulging scenarios\n(a>3*s/4) exhibit a negative peak first and then turn to a positive value when it comes to the\ncomplete collision. However, we would like to point out that the displacement of the mass\ncenter in east-bulging models is relatively small. The original calculation with a fixed mass\ncenter should make more sense. The duration of a negative torque and its magnitude should\nbe minor (Fig. R3c,d). Besides, the overriding Eurasian plate is fixed on the west side and free\non the east side in the models. The asymmetry-induced clockwise rotation of the Eurasian\nplate could facilitate the accumulation of the eastern oceanic slab (Extended Data Fig. 7), and\nthe resulting greater slab pull would weaken the initial clockwise rotation. Therefore, the\nclockwise rotation is not significant in the east-bulging numerical model. And the first-order\nphenomenon, the absence of positive peaks of torques for the east-bulging scenarios (a>s/2),\nis consistent with the disappearance of the first peak in rotation rates, which precludes the\npossibility that the first collision occurred on the eastern side.\nFig. R3. Calculation of the resistive torques. a, Diagram of the Indian continent before collision.\nThe center of mass of the Indian plate is set to move westward during the collision process.\nOther captions are the same as in Extended Data Figure 5. b, Calculated resistive torque versus\ncontact width between continents. c and d are corresponding diagrams and curves considering\nthe center of mass remains fixed.\nFurther comments:\nline 114: But, North-South convergence, as described in the Methods (l. 499/500) is not\nidentical to E-W boundary, unless the convergence is orthogonal to the boundary. Is that the\ncase?\nA short answer is Yes. We have considered this when constructing the model. Here, we extracted\nthe azimuth of the Cenozoic trajectory of the central reference point, which represents the\ndrifting direction of the Indian plate. The average of azimuth spans from 40° to 20° (Fig. R4),\nwhich is consistent with that derived by Gibbons et al. (2015). We have also calculated the\nstrike of the trench derived from different proxies. The recent paleomagnetism studies\nconstrains that the azimuth spans from 317° to 290° (Yi et al., 2015; Tong et al., 2022; Bian\net al., 2022; Ma et al., 2022). The restoration of the intra-Asian shortening and extrusion\nshows that the boundary is oriented at around 293°-291° (van Hinsbergen et al., 2019).\nThe azimuth of the linear high-velocity seismic anomaly in the lower mantle is around 301 ° -305\n° (Replumaz et al., 2014). Overall, the convergence boundary before the collision is oriented\nat around 290°-310°. Thus, the convergence is almost orthogonal to the boundary (Fig. R4).\nFor clarification, we have added extra descriptions in the Methods (lines 117-121, 551-554).\nFig. R4. Convergence direction of the Indian plate derived from the four reconstructions (left)\nand the diagram of the convergence direction and the azimuth of the convergence boundary\n(right).\nline 238: Isn't the pre-dating by 5-10 Myrs best seen in Fig. 3a?\nAgree. We have added Fig. 3a as an additional figure reference. See line 244.\nline 265: Why is Sibumasu Block abbreviated with \"SP\"? What does the \"P\" stand for?\nThe \"SP\" here represents Shan Plateau, another widely-used name for this terrane. We have\nchanged the abbreviation here to \"SB\". See line 271.\nline 482: \"as long as the Euler Poles do not overlap with the North Pole\" - I am not sure what\nyou mean with this. I mean, also if there is an Euler pole at the equator 90° away from the\ncontinent, this corresponds to pure northward motion with no rotation.\nThanks for your comment. We have corrected it here. See lines 532-533.\nFor the absolute viscosity values in line 508, do you mean in reality or in your model?\nThe values here refer to the absolute viscosity in our models. We have added “in the models”\nin this sentence. See line 562.\nI am also wondering whether you need to specify a higher viscosity in the lithosphere; doesn't\nit come in automatically with the higher (lower?) temperature, so you don't need to explicitly\nprescribe a higher viscosity (because otherwise, the same effect is considered twice)?\nThanks for pointing this out. We have added extra explanations in the Methods (lines 570-\n572). The viscosity of the oceanic mantle lithosphere is generally controlled by temperature in\nthe models, as its temperature and composition (olivine) are relatively well known. But for the\ncontinent, its geotherm and compositions are more complex. For simplicity, we mainly\nconsider the effective viscosity for continents, whose value could be easily approximated by\nthe composition tracers of the crust and lithosphere mantle. Therefore, we use a special\nviscosity setting in the continent, which means that if an element contains composition tracers\nof continental lithosphere, then its composition-dependent viscosity will take over and the\ntemperature effect will not take effect. In theory, the temperature dependence should also\naffect viscosity, but this is neglected for the above purpose.\nIn your viscosity law in line 512, everything seems to be dimensionless, so how about the\nsecond invariant of the strain rate tensor, this should normally be in units of 1/s; how do you\nget this dimensionless?\nYes, the viscosity law here is dimensionless. In the calculation of the second invariant of the\nstrain rate, every parameter including velocity and length is dimensionless, which makes the\nsecond invariant of the strain rate dimensionless automatically. The dimensionless factor of\ntime equals 𝑅0\"/𝜅0, in which 𝑅\n0\nis the radius of the Earth, and 𝜅\n0\nis the reference thermal\ndiffusivity. But if you consider the dimension in the viscosity law, the unit 1/s of the strain\nrate tensor is canceled out by the unit 𝑠(1_n)tn in the prefactor A(r).\nline 571: I think the sentence \"The width ...\" can be improved. I prefer how it is said on line\n109/110 \"oceanic lithosphere on the eastern side of the Indian continent is much wider than\nthat on the west\"\nWe have changed the expressions here. See lines 648-650.\nExtended Data Figure 5: Why is there a \"1e16\" on top of panel b? Why you don't include\ncurves for a=0 and a=1?\nThe 1e16 here is the multiplier of the numbers beside the vertical axis. We have redrawn\nExtended Data Figure 5 and included the curves for a=0 and a=1.\nExtended Data Table 3: So here you don't distinguish between Upper Crust, Lower Crust and\nLithospheric Mantle? And the distinction between Asian Non-Cratonic and Cratonic\nLithoshere is for all three (Upper Crust, Lower Crust and Lithospheric Mantle)?\nWe have considered the viscosity layering of the continental lithosphere in the original models\nand corrected the chart accordingly (Extended Data Table 3). The distinction of viscosity\nbetween Asian non-cratonic and cratonic lithosphere is for the whole lithosphere.\nExtended Data Table 4: Does \"Neutral\" mean it is the same as the underlying mantle, apart\nfrom Temperature dependence? And what is actually the temperature dependence of density\n(i.e. thermal expansivity)? I don't see it anywhere specified.\nWe have added the description of density in the Methods (lines 596-603). Temperature-,\ncomposition-, and phase-change-dependent density is considered in the numerical modeling.\nThe word \"Neutral\" here means that the composition-dependent density of the lithospheric\nmantle is the same as that of the underlying mantle. The thermal density anomaly is defined as\n-a 0 , 0 LT. The a 0 in this expression is the thermal expansivity, which equals 3.0*10-5 /K in\nthe modeling.\nMinor comment:\nline 33: \"has been debated\"\nThere are two subjects in the sentence, i.e. \"when India collided with Eurasia\" and \"how this\nprocess varied along the strike of the Indian front\". Thus we think the word \"have\" here should\nbe more appropriate.\nline 294: \"could have either\"\nThank you. This has been changed. See lines 300-301.\nReviewer #2 (Remarks to the Author):\nThis is an interesting paper bringing some new constraints to understanding the history and\npotential driving forces of the India-Eurasia collision. Using the connection between changes\nin the azimuth of the India-Eurasia collision to constrain geodynamic models is a novel way of\nlooking at this problem. I enjoyed reading the paper, and suggest that it should be published\nwith some relatively minor corrections. The English language needs a bit of improvement\nthroughout the paper. I have made some comments along these lines.\nThank you very much for your positive comments and valuable suggestions on our paper that\nhave greatly improved our manuscript. We have also improved the English language based on\nyour comments.\nIn terms of the interpretation of the results, I have this overall comment: The authors are\nfollowing the premise that the time-dependent changes in rotation of India relative to Eurasia\ncan only be explained by a successive continent-continent collision which needs to start around\n55Ma or earlier, necessitating an at least 2000 km wide thinned continental margin off\nnorthwestern India (the details of the inferred geometry are difficult to extract from the paper).\nWe appreciate this brief summary of the reviewer. This is a good point that necessitates a brief\ndiscussion of the different modes of India-Eurasia collision, as well as the morphology and\nposition of Greater India during collision, as part of the broad implication of this study.\nInitially, we mainly focus on the NW-SE asymmetry-induced time-varying resistive torque\nand rotation rate. Then we further discuss different collision scenarios, where we show that\nthe NE-SW width, crustal thickness, and characteristics (continental or oceanic) of Greater\nIndia are subject to heavy debate and that our model is compatible with multiple hypotheses.\n(1) The first question one needs to ask oneself is: Is this geodynamically reasonable or possible?\nThe widest continental margin in existence today is about 400 km in width. There is a simple\nreason for this. It is not rheologically possible to create a much wider passive margin, as the\nstretched lithosphere will simply break, producing seafloor spreading, preventing excessive\nstretching from occurring. Stretching factors larger than four lead to SFS (see classical papers by\nDan McKenzie). Continental crust is too weak to be stretched beyond this limit. So how exactly\ndo you propose to generate a 2000 km wide passive continental margin to start with, 5 times\nwider than any margin in existence? In my view this essentially falls in the land of geodynamic\nphantasy, so I'd like to see this justified in the paper (personally, I don't think this is possible, but\nI am looking forward to see the explanation you come up with!).\nThis is a good question. We agree with the reviewer that geodynamically it is difficult to\ngenerate 2000-km wide thinned continental crust through a single stretching event. However,\nsuch a continental crust could likely be generated via multiple instead of a single extension\nevent, during which the underlying mantle lithosphere could temporally restabilize/regrow\n(McKenzie, 1978) to prevent necking and detachment as occurs in a single stretch. Since the\ncrustal volume remains largely unchanged, its thickness may keep decreasing and length\nincreasing. These extension events might have occurred during the ~280 Ma eruption of Panjal\nplume in northwestern Greater India and southern Qiangtang (Shellnutt et al., 2011; Zhai et al.,\n2013), the ~220 Ma breakup of Lhasa terrane next to northeastern Greater India and\nnorthwestern Australia caused by southward subduction (Zhu et al., 2011; Xu et al., 2022),\nthe ~155 Ma eruption of plume in northeastern Greater India and northwestern Australia\n(Ludden, 1992; Wang et al., 2020), the ~130 Ma eruption of Comei plume in northeastern Greater India\n(Zhu et al., 2009), and the significant slab pull caused by double subduction during late\nCretaceous (Jagoutz et al., 2015). The widely distributed extension events not only provided\nvertical plume push or lateral tensile force (Brune et al., 2023), but also heated and weakened\nthe original lithosphere, which facilitate the formation of a wide stretched lithosphere (van der\nPluijm & Marshak, 2004).\nOn the other hand, our opinion remains neutral in terms of the crustal thickness of the Greater\nIndia. From our understanding, the continental crust with a normal or thinner thickness could\nboth be possible scenarios. The former scenario has been widely accepted by previous studies\n(e.g. Ingalls et al., 2016; Ding et al., 2022). The latter one has been proposed by several recent\nstudies (Himalandia in Liu et al. (2023); Zealandia-type in Li & Robinson, (2023)) through\nnumerical modeling and mass-balance calculation. They think Greater India was similar to\neastern Asia and southwestern Pacific nowadays, with a wide (>1000 km) and thinned\ncontinental crust as revealed by seismic probing (Gallais et al., 2019). Besides, we do not\npreclude the scenario that there were some seafloors spanning across Greater India induced by\nstretching as suggested by van Hinsbergen et al. (2017) and Yuan et al. (2021). We suggest that\nthese would not change the major rotation features, since the stretched continental lithosphere and\nthe new oceanic lithosphere are still more buoyant than the surrounding old oceanic lithosphere\n(Liu et al., 2023), thus could provide continuous resistive forces and torques altering rotation\nrates, similar to those discussed in the main text. In this paper, we remain neutral to these\nproposals and we have added related discussions in the main text (lines 351-357).\nThe subduction of this much continental crust also presents other problems, as is discussed,\nfor instance, in van Hinsbergen DJ, Lippert PC, and Huang W. Unfeasible subduction? Nature\nGeoscience. 2017 Dec;10:878-9.\nWe have carefully looked through the related papers by van Hinsbergen, Rowley, and their\ncolleagues (Ingalls et al., 2016; van Hinsbergen et al., 2017; Rowley & Ingalls, 2017). While\nvan Hinsbergen et al. propose that there is a mismatch between the kinematic convergence and\ntectonic recovery and only the \"Greater Indian Basin\" model could account for the discrepancy\n(van Hinsbergen et al., 2017), Rowley et al. insist that a large amount of erosion of the upper\ncrust should be considered and thus could fill the gap, and that the classic \"Greater India\"\nmodel still works (Ingalls et al., 2016; Rowley & Ingalls, 2017). Besides seafloor spreading and\nsignificant erosion, there are several plausible explanations for the mass discrepancy of the\ncontinental crust. First, a portion of the felsic upper crust of the normal Greater India might be\ndragged by the preceding slab and underwent phase change-induced densification. The\nnegative buoyancy would make them subduct to the deeper mantle instead of exhumed to the\nground (Wang et al., 2022). Second, there could also be a Greater India with a thinner\ncontinental crust as illustrated above. Thus less amount of the upper crust might be scrapped\noff and preserved nowadays (Liu et al., 2023; Li & Robinson, 2023). As discussed above,\nGreater India with a normal or stretched lithosphere are both plausible scenarios in fitting\nrotation features. In this paper, we remain neutral to these ideas and we have added related\ndiscussions in the main text (lines 351-357).\n(2) Along the same lines, I am surprised that the recent paper by Pusok and Stegman (The\nconvergence history of India-Eurasia records multiple subduction dynamics processes, Sci\nAdv 2020) is not cited at all. Essentially it provides some ideas for an alternative way to create\na two-stage collision history, with the first collision reflecting the collision of an island arc\nwith Eurasia followed later by a continent-continent collision. This scenario was also\nproposed in a number of other papers. This should at least be discussed as a potential\nalternative way to explain the kinematic observations.\nThank you for the reminder. We have added related discussions in the main text (lines 337-\n351). We have read related papers and here is our analysis. In short, the earlier \"Arc-Continent\nCollision\" scenario is compatible with the current collision process. According to the regional\nplate reconstruction (Zahirovic et al., 2016), the eastern portion of the India-Eurasia\nconvergence boundary was a single subduction zone after the Woyla arc collided with\nSumatra at ca. 95 Ma (Advokaat et al., 2018). The double subduction zone is restricted to the\nwestern portion, consisting of Ladakh-Kohistan arc in the west and West Burma terrane in the\neast (Westerweel et al., 2019; Martin et al., 2020). It has been proposed that double\nsubduction could exert a larger force on the subducting plate and thus could account for the\nhigh drifting rate of the Indian plate during the Cretaceous-Paleocene boundary (Jagoutz et al.,\n2015; Pusok et al., 2020). Considering the width of the western and the eastern convergence\nboundary was almost the same, there would be a larger slab pull in the west than that in the\neast (Fig. R5a), which would generate prominent torque resulting in extra clockwise rotation\nof the Indian plate. The torque and the associated extra clockwise rotation would come to an\nend when the island arc collided with either Indian or Eurasian continent, resulting in the\nextinction of the asymmetric slab pull along the boundary (Fig. R5b).\n[Redacted].\nBut the time of the arc-continent collision is still controversial (Najman et al., 2017). Several\nstudies propose that the arc-continent collision had been finished before 60 Ma, such as 10080\nMa (Zaman et al., 2013), 90-80 Ma (Rehman et al., 2011), 85 Ma (Chatterjee et al., 2013;\nBorneman et al., 2015), and 70-60 Ma (Burg et al., 2011). If this is the case, the arc-continent\ncollision and its associated cessation of extra clockwise plate rotation should be the overture\nof the Paleocene continent-continent collision. The rotation rate during the late Cretaceous,\nhowever, is heavily influenced by plume push (Extended Data Fig. 8). Thus it is difficult to\nisolate the signal of arc-continent collision directly.\nHowever, some studies believe that the arc-continent collision happened at ca. 60-50 Ma\n(Aitchison et al., 2007; Khan et al., 2009; Hébert et al., 2012; Bouilhol et al., 2013), which are\nfurther cited in some recent paleomagnetism and numerical modeling studies (Jagoutz et al.,\n2015; Westerweel et al., 2019; Martin et al., 2020; Pusok et al., 2020). If this is true, one may\nthink that the rotation peak at 52-44 Ma might correspond to the arc-continent collision, and the\npeak at 33-20 Ma might reflect the continent-continent collision. But the earlier cessation of the\nextra clockwise rotation and the later onset of the extra counterclockwise rotation should\ndisplay a continuously increasing rotation rate curve from ca. 60-50Ma to ca. 30-20 Ma (Fig.\nR6a), which is not consistent with the Cenozoic double peak observation (Figs. 1c, R6b).\nOverall, we think the arc-continent collision before 60 Ma could be a possible scenario of the\nearlier collision process, which is compatible with the current Paleocene collision process.\nFig. R6. Analogue extra rotation rate of the two collision scenarios. a, Arc-continent collision\nat ca. 55 Ma. b, Greater India collision at ca. 55 Ma proposed in the main text.\nSome detailed comments:\nLine 16:\nThe Cenozoic India-Eurasia collision has profound impacts on the formation of the\nmodern Tibetan plateau.\nshould be\nThe Cenozoic India-Eurasia collision has had profound impacts on the formation of the\nmodern Tibetan plateau.\nChanged. See line 16.\nLine 136, and many subsequent instances:\nI do not understand why in \" West-bulging model\", \"Middle-bulging model\" etc \"west\" and\n\"middle\" is capitalized. It should be lower-case.\nWe have changed the expressions throughout the main text.\nLine 145:\nThe two modeled rotation peaks postdate the respective initial collision and complete\ncollision by about 8 Myrs (Fig. 2a-d), which bears their formation mechanism.\nshould be\nThe two modeled rotation peaks postdate the respective initial collision and complete\ncollision by about 8 Myrs (Fig. 2a-d), reflecting their formation mechanism.\nChanged. See line 150.\nLine 150: neglectable\nshould be\nnegligible\nChanged. See line 154.\nLine 154: This results in no noticeable rotation.\nThe use of the word rotation is ambiguous here, as well as through much of the rest of the\npaper. The issue is that all plates rotate about Euler poles at all times, ie plates are always in\nrotation. What you mean here is rotating about a pole located proximal to the Indian continent.\nYou need to think about how to do that better, ie you cannot write that at some point India was\nnot rotation, and then it started to rotate, seeing that it is ALWAYS rotating.\nThanks for the suggestion. This is a good point. We actually have also realized this. \"rotation\"\nis a more convenient word in many different circumstances throughout the manuscript, but the\nactual meaning might be better expressed as “azimuthal change” of a plate. We now add this\ndefinition in the abstract (line 20). And we have added the word \"extra\" before \"rotation\" in\nline 159 and line 203.\nLine 166: Similarly, a \"peak in rotation\" is not the correct term here. You actually mean a\npeak in rotation rate. The same issue occurs in other places in the text.\nThanks again for the suggestions. We have changed the expressions throughout the main text.\nLine 261: Some of the letter labels overlying \"Greater India\" are barely visible, especially\nthose frequently referred to in the text like the \"Saga\" and \"Tingri\" areas.\nWe have redrawn Figure 4 to make it clearer."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have done quite an effort to revise the paper, and all my comments have been addressed. Accordingly, I have only a few small comments - partly things I had overlooked in the first round, partly new issues that came with the revisions.\nOn lines 94-98 you write now that the first peak is not as robust and that this is likely due to the discrepancy in the adopted absolute reference frames. I suggest that, like you explain in the rebuttal letter, you can, for \"evidence\", add a sentence that, if you replace the slab-fitting reference frame by a hotspot reference frame, you still get a 52-44 Ma peak.\nAnother issue is about the shift in the center of mass, which is newly introduced. Is this just something that you prescribe (and which is the same in Figures 2 and Extended Data Figures 5a and 6, and not considered in Extended data Figures 5c and 7) or do you actually compute it for all your models? You explain why it should move westward, but shouldn't it also move southward relative to the continent, because material is lost in the north? Or do you only consider the westward component of motion? These questions apply to both your numerical models and your simplified quantitative calculation.\nRelated to this, you refer to Extended Data Figure 5d on line 205. You explain in your rebuttal letter and the methods section that the east-bulging case is better combined with no motion of the center of mass, but here it does not become clear why you refer to panel d rather than b.\nCaption Figure 3: I think you should add to the figure caption, like you did in your rebuttal letter, that the dashed lines outline the subducted portion.\nline 345: You mean, a continuously increasing counterclockwise rotation rate, i.e. that the rotation rate was always counterclockwise, despite the additional clockwise rotation?\nlines 531-533: I would rewrite this sentence as follows: \"... which has an inherent azimuth variation component, except in special circumstances, such as when the Euler pole is located at the North/South pole, or at the equator 90° of longitude away from the central reference point\"\nlines 570-572: You write in your rebuttal letter that temperature is not considered in the continental lithosphere, but you don't write it explicitly here. If that is so, I think you should also write it here for clarity.\nExtended data Table 3: You added eta_cut-off, but I think you need some explanation that these are the upper and lower limits of viscosity considered. Is viscosity given in units of 10**21 Pas? In general, I find the mix of dimensionless and dimensional numbers a bit confusing. You need to be clear about what is what - for dimensional numbers give the units, and for dimensionless ones give the normalizing (or reference) value.\nMinor comments:\nline 146: perhaps better \"timing\" rather than \"tempo\"?\nline 204: add \"simplified\" before \"quantitative calculation\"?\nline 220: better \"abnormal\"?\nline 635: remove \"peak\"?\nline 639: instead of \"is larger\" perhaps better \"reaches a larger value\", since it is not always larger?", "reviewer_2": "The authors have done an excellent job in revising the paper, addressing all reviewer comments in detail. I am happy for the paper to be published in its present form."}, "author_response": "Reviewer #1 (Remarks to the Author):\nThe authors have done quite an effort to revise the paper, and all my comments have been\naddressed. Accordingly, I have only a few small comments - partly things I had overlooked in\nthe first round, partly new issues that came with the revisions.\nThank you very much for the supportive comments and valuable suggestions that have greatly\nimproved our manuscript. The remaining issues have been corrected now.\nOn lines 94-98 you write now that the first peak is not as robust and that this is likely due to the\ndiscrepancy in the adopted absolute reference frames. I suggest that, like you explain in the\nrebuttal letter, you can, for \"evidence\", add a sentence that, if you replace the slab-fitting\nreference frame by a hotspot reference frame, you still get a 52-44 Ma peak.\nThank you for your suggestion. We have added the related sentence in the main text. See lines\n102-104 in inline markup mode.\nAnother issue is about the shift in the center of mass, which is newly introduced. Is this just\nsomething that you prescribe (and which is the same in Figures 2 and Extended Data Figures\n5a and 6, and not considered in Extended data Figures 5c and 7) or do you actually compute it\nfor all your models? You explain why it should move westward, but shouldn't it also move\nsouthward relative to the continent, because material is lost in the north? Or do you only\nconsider the westward component of motion? These questions apply to both your numerical\nmodels and your simplified quantitative calculation.\nThanks for pointing this out. Since the exact movement of mass center in the models is difficult\nto extract, what’s presented is a simplified version of the movement of mass center\ncorresponding to the simplified torque analysis. As the numerical models are designed to be\nfree-slip and without any prescribed forces on the boundaries, the rotation of the subducting\nplate and the shift of the center of mass are both natural outcomes of the model without\nadditional prescription. In order to explain the observed rotation feature of the numerical\nmodels, we introduced the movement of the center of mass in the simplified torque calculations,\na phenomenon which is observed in the west- and middle-bulging models accompanying\nsignificant counterclockwise rotations. As for the east-bulging model, the magnitude of\ncounterclockwise rotation and the offset of the center of mass is relatively minor, therefore the\ntorque calculations with a fixed mass center should be more appropriate.\nWe agree with the reviewer that the center of mass should move southward with the collision\ngoing on. We have redrawn the figures to fix this problem (Fig. 2, Supplementary Figs. 6,7).\nActually, these changes do not influence the original explanation because the length of the lever\narm is unaltered (Fig. 2). Therefore the southward movement of the center of mass does not\nchange the resistive torques in the simplified calculations. We have added a related explanation\nin Methods. See lines 474-476.\nRelated to this, you refer to Extended Data Figure 5d on line 205. You explain in your rebuttal\nletter and the methods section that the east-bulging case is better combined with no motion of\nthe center of mass, but here it does not become clear why you refer to panel d rather than b.\nThanks for pointing this out. We have made some clarifications in the figure caption that can\nguide the readers to understand the different choices of torque analysis for different models by\nreferring to the Methods part. See figure captions below Supplementary Figure 5 in the separate\nsupplementary file. We avoid to discuss the movement of mass center in the main text to keep\nthe main text concise and to avoid redundancy, since all related discussions are compactly\nshown in the Methods section.\nCaption Figure 3: I think you should add to the figure caption, like you did in your rebuttal\nletter, that the dashed lines outline the subducted portion.\nThank you for your suggestion. We have added the related sentence in the caption of Figure 3.\nSee lines 222-223.\nline 345: You mean, a continuously increasing counterclockwise rotation rate, i.e. that the\nrotation rate was always counterclockwise, despite the additional clockwise rotation?\nHere we use the expression \"a continuously increasing rotation rate\" to describe the\nhypothetical variation in the rotation rate of the Indian plate, i.e. the earlier weakening of the\nclockwise rotation induced by arc-continent collision, and the latter strengthening of\ncounterclockwise rotation induced by continent-continent collision. There could be a period of\nnegative (i.e. clockwise) rotation rate when there was a double-subduction system exerting\nlarge enough extra clockwise torque on the Indian plate. But the observed Cenozoic rotation\nrate of the Indian plate was positive (i.e. counterclockwise) almost all the time (Fig. 1c,\nSupplementary Fig. 2), and the double-peak feature also seems inconsistent with the predictions.\nTherefore the scenario in which the arc-continent happened at ca. 60-50 Ma is not favored.\nWe now change this expression to \"a continuously increasing rotation rate (assuming negative\nfor clockwise rotation and positive for counterclockwise rotation)\". See lines 353-354.\nlines 531-533: I would rewrite this sentence as follows: \"... which has an inherent azimuth\nvariation component, except in special circumstances, such as when the Euler pole is located at\nthe North/South pole, or at the equator 90° of longitude away from the central reference point\"\nThank you for your suggestion. We have changed the expression here. See lines 387-389.\nlines 570-572: You write in your rebuttal letter that temperature is not considered in the\ncontinental lithosphere, but you don't write it explicitly here. If that is so, I think you should\nalso write it here for clarity.\nThank you for your suggestion. We have added related description in the Methods. See line 429.\nExtended data Table 3: You added eta_cut-off, but I think you need some explanation that these\nare the upper and lower limits of viscosity considered. Is viscosity given in units of 10**21 Pas?\nIn general, I find the mix of dimensionless and dimensional numbers a bit confusing. You need\nto be clear about what is what - for dimensional numbers give the units, and for dimensionless\nones give the normalizing (or reference) value.\nThank you for your suggestion. We have added related descriptions below the Supplementary\nTable 3. See table caption below the Supplementary Table 3. We have also converted all non-\ndimensional values to dimensional values, including those in Supplementary Tables 2 and 3.\nMinor comments:\nline 146: perhaps better \"timing\" rather than \"tempo\"?\nline 204: add \"simplified\" before \"quantitative calculation\"?\nline 220: better \"abnormal\"?\nline 635: remove \"peak\"?\nline 639: instead of \"is larger\" perhaps better \"reaches a larger value\", since it is not always\nlarger?\nThank you again for your suggestions. We have changed related expressions in the main text\nand Methods. See lines 152, 209, 227, 494, 499.\nReviewer #2 (Remarks to the Author):\nThe authors have done an excellent job in revising the paper, addressing all reviewer comments\nin detail. I am happy for the paper to be published in its present form.\nWe appreciate your thorough review, which has greatly improved our manuscript."}]

Decision letter