Cenozoic India-Asia collision driven by mantle dragging the cratonic root

Domain: earth_science
Published: 2024-08-06
Editorial outcome: major_revision

Abstract
The driving force behind the Cenozoic India-Asia collision remains elusive. Using global-scale geodynamic modeling, we find that the continuous motion of the Indian plate is driven by a prominent upper-mantle flow pushing the thick Indian lithospheric root, originated from the northward rollover of the detached Neo-Tethyan slab and sinking slabs below East Asia. The maximum mantle drag occurs within the strong Indian lithosphere and is comparable in magnitude to that of slab pull (10^13 Nm^-1). The thick cratonic root enhances both lithosphere-asthenosphere coupling and upper-plate compressional stress, thereby sustaining the topography of Tibetan Plateau. We show that the calculated resistant force from the India-Asia plate boundary is also close to that due to the gravitational potential energy of Tibetan Plateau. Here, we demonstrate that this mantle flow is key for the formation of the Tibetan Plateau and represents part of a hemispheric convergent flow pattern centered on central Asia.

Introduction
The continuous northward convergence of the Indian plate with Tibet after their Cenozoic (~50 Ma) collision represents a key linkage among the termination of the Neo-Tethyan ocean, the birth of the Indian Ocean, and the formation of Tibetan Plateau. In addition, geophysical analysis suggests that the late-Cenozoic resistance force due to the gravitational potential energy (GPE) of the Himalayas and Tibetan Plateau ranges from 3×10^12 Nm^-1 to 8×10^12 Nm^-1. Consequently, the force driving the northward convergence should be only higher than this value given the potential resistance from other boundaries of the Indian plate. Yet, direct slab pull, the most commonly invoked plate driving force with a comparable amplitude (10^13 Nm^-1), is largely negligible because of the absence of a continuous oceanic slab below the present collision zone. Ridge push, on the other hand, is too small, given its typical amplitude of 3×10^12 Nm^-1. Thus, the post-collisional Indian motion cannot be adequately explained by the above mentioned plate driving forces, leading to the intriguing question of what the dominant driving mechanism and its associated origin are.

This problem has invoked multiple hypotheses. One of these refers to the continuous subduction of the continental lithosphere which is considered dense enough to provide an alternative slab pull force. However, their predicted depth-continuous steep slab is at odds with the observed seismic structure where the Indian lithosphere is shallowly underplating Eurasia with no connection to prominent deeper slabs, not to mention the likely early formation of eclogite densification and the early slab breakoff. Another proposal suggests that nearby subducting plates like the Capricorn and Australian plates transmit far-field slab pull from below Southeast Asia to the Indian plate interior, but this mechanism is dynamically challenging for neighboring plates with near-parallel motion directions, especially for a force large enough to overcome the substantial GPE-induced resistance from Tibet.

A conceptually more feasible driving force is the mantle drag from the faster-moving asthenosphere acting at the bottom of the lithosphere, where the associated fast mantle flow was proposed to be induced either by a rising plume and/or suction of foundered slabs. However, many aspects of this driving force remain unclear, such as its amplitude and duration, how the asthenosphere interacts with the lithosphere, the physical origin of the fast mantle flow, and its capability to maintain both plate motion and the India–Asia collision.

In this work, we address these questions using a global, time-dependent geodynamic model that couples the lithosphere and the convective mantle. We further quantify both the driving and resisting forces within the model. Through a systematic and quantitative analysis of the entire dynamics system, we demonstrate that mantle drag is the dominant driving force behind the Cenozoic India–Asia collision and the motion of the Indian plate.

Results
Cenozoic slab evolution and mantle flow
We tackle these problems using a global, sequential-in-time data-assimilation model starting from 200 Ma (see Method for more details) based on a recent plate reconstruction. This model generates evolving mantle structures (Supplementary Fig. 1) including slabs and large-low-shear-velocity provinces (LLSVPs) that match present-day tomographic structures both regionally (see below) and globally. Its temporally evolving continental lithospheric thickness in response to mantle dynamics also allows quantifying the lithosphere-asthenosphere interaction due to the variable lithosphere-asthenosphere boundary (LAB) depth (Supplementary Fig. 1c, d), as recently done. To reproduce observed plate motion at different times using a fully dynamic approach, we solve the flow using the corresponding mantle structures from the data-assimilation model while allowing the surface to deform freely in the horizontal directions.

During this calculation, we also vary the effective viscosity values of the plate boundaries and the oceanic asthenosphere, both of which are not well determined during the initial sequential simulation. After obtaining a satisfactory match between modeled and observed plate motion at the present, the same viscosity adjustment is applied to all other times (see Method). Finally, based on these results we analyze the properties and origin of the relevant forces.

Our model’s reproduction of present-day global plate motion (Fig. 2a) is comparable to, or better than published results. As a step further, our calculation also properly reproduces the past motions of the Indian plate (Figs. 2b, c, and 4c). A prominent dynamic feature is that the asthenospheric flow south of the Indian continent is consistently much faster than the surface plate motion. The good fit between modeled (see Method) and observed asthenospheric (200 km depth) seismic anisotropy (Fig. 2d) further supports the directionality of this fast mantle flow. These results indicate that the previously hard-to-quantify viscous mantle drag should act as an important driving force for the persistent Cenozoic Indian plate motion.

In our model, the faster-than-surface asthenospheric flow is a direct result of the lateral dynamic pressure gradient in the upper mantle, where the pressure is high beneath the southern Indian Ocean and low beneath the Indian continent (Fig. 3b). This pressure gradient exerts a net force pushing the mantle to go northward, generating a type of flow called the Poiseuille flow. Consequently, the overriding Indian lithosphere, which moves more slowly due to resistance from its boundaries, especially on the northern side, experiences a mantle drag from the fast-moving asthenosphere below it.

The low-pressure upper mantle beneath the Indian continent is a combined result of (1) the sinking and southward rollover of the Neo-Tethyan slab in the mid-lower mantle, and (2) the voluminous slab pile beneath Eurasia at lowermost mantle depths comprising both the Neo-Tethyan and Izanagi slabs. The southward rollover of the sinking Neo-Tethyan slab at lower-mantle depth, caused by continuous trench advance, reduces the pressure beneath the subducting plate instead of the overriding plate, with the latter being more commonly observed with retreating slabs. The relatively high pressure in the upper mantle south of the continent is partly due to the upwelling plume from the underlying African LLSVP. Both the subducted slabs and upwelling plumes are consistent with seismic tomography (Fig. 3a).

The amplitude of dynamic pressure and its lateral gradient scale linearly with the density anomalies and only weakly depend on mantle viscosity, for a given spatial pattern of mass distribution. This property allows us to estimate the relative contributions of various density anomalies to the observed present-day plate motion and associated mantle drag, without worrying about the effect of mantle viscosity. Compared to the Eurasian slab pile, the Neo-Tethyan slab is less voluminous, thus its removal only slows the Indian plate by ~20% (Supplementary Figs. 2c, 3b). However, in this case, the Eurasian plate moves northward much faster at a similar rate to India, implying little convergence between the two plates. Removing lower mantle slabs globally will stop the mantle flow and its driving of the Indian plate (Supplementary Figs. 2f, 3c), consistent with the scenario where both slab pull and slab suction are absent. Thus, both slab piles are essential in driving Indian convergence: while the voluminous Eurasian slab pile generates regional-scale low pressure below central-eastern Asia, the Neo-Tethyan slab interjects to expand the low-pressure center to below north India and generate the proper asthenosphere flow pattern associated with the observed India–Tibet convergence (Fig. 3b). The dense continental eclogite (see below) existing in all three cases have only minor impact on the large-scale mantle flow. This highlights a complex multiscale interaction among the sinking mantle slabs that is hard to illustrate in simpler models used previously.

As a further step from earlier studies, the time-dependent nature of our model allows us to examine the evolution of the slab and its control on the mantle flow during the geological past. The detachment of the Neo-Tethyan slab happened between 50 Ma and 40 Ma after the closure of the Neo-Tethyan ocean (Fig. 3d, e). In contrast to the traditional slab suction hypothesis, the faster-than-surface upper-mantle flow appeared even before the slab detachment (Fig. 3d). This flow was likely attracted by the Izanagi subduction in the north, which also induced the Neo-Tethyan trench advance and slab rollover. This fast asthenospheric flow persists from >50 Ma to the present day, suggesting a consistent driving force for the Indian plate throughout the post-collision Cenozoic era. In the following section, we quantify this driving force and demonstrate its dominant contribution to the India–Asia collision.

Mantle drag due to lithosphere-asthenosphere interaction
While some previous research suggested mantle flow as a potential driving mechanism for the post-collisional Indian motion, the existence, magnitude, distribution, and evolution of this force all remain unclear. For example, the much lower (by many orders of magnitude) asthenospheric viscosity than that of the overlying lithosphere seems to dismiss mantle drag as a major driving force. Thicker lithospheric keels could affect the coupling between the lithosphere and the underlying mantle, but this effect on continental motion is uncertain, due to its strong dependence on the poorly constrained mantle flow pattern. Therefore, it is necessary to further quantify the flow-induced mantle drag and how the asthenosphere interacts with the lithosphere with the presence of a weak asthenosphere and variable LAB depth in an evolutionary model.

Mantle drag consists of two types (Fig. 4a, b): (1) Friction drag due to shear stress along the lithosphere boundary, and (2) pressure drag due to force normal to the contact surface through pressure gradients. Inside the convecting Earth, friction drag originates from the vertical gradient of horizontal flow around the base of the lithosphere, while pressure drag mostly occurs along lithospheric steps where the LAB depth varies laterally. Our analysis (Method) shows that the magnitude of total mantle drag along the plate motion direction is 1.1×10^13 Nm^-1 at present, exceeding the highest estimation (8×10^12 Nm^-1) of Tibetan resistance from the northern plate boundary. Among these forces, friction drag has a dominant magnitude of 8.8×10^12 Nm^-1 while pressure drag is 2.2×10^12 Nm^-1 (Fig. 4a, b).

Interestingly, both friction drag and pressure drag show a prominent concentration around the southwestern edge of the Indian continent, where a large lithospheric step exists. This concentration accounts for ~50% of the total friction drag and nearly 100% of the pressure drag. Spatially, this lithospheric step corresponds to a sharp slowdown of the asthenospheric flow that is blocked by the thicker continental keel to the northeast (Fig. 2b). This observation confirms the important role of a thicker lithospheric root for increasing coupling between the plate and the asthenosphere.

The resistant force acting on the northern plate boundary is 9.7×10^12 Nm^-1 at present-day (Fig. 4b), comparable to the total mantle drag as well as the GPE estimation (3×10^12 Nm^-1 to 8×10^12 Nm^-1). This means that most of the driving force is balanced by the collisional resistance from the India–Asia collision. The largest resistant force is located in the mid-western part of the collision zone, likely as a direct reaction along the plate motion direction to the large driving force from the southwestern edge of the Indian continent. This shows that our model, although not explicitly including a free surface to mimic the high plateau topography, captures the correct sub-surface vertical stress along the plate boundary that represents the correct GPE as occurred in reality.

In contrast to the traditional view that a weak asthenosphere diminishes the mantle drag, our results demonstrate that the maximum shear stress occurs within the relatively strong lithospheric root (Supplementary Fig. 4a). Inside the weak asthenosphere, the maximum friction drag is about half of that within the lithosphere above (Supplementary Fig. 4b, c), much greater than that implied from a linear proportion scaled with viscosity, as previously found. In general, the effective mantle drag gradually increases from the weak asthenosphere to the strong lithospheric interior, irrespective of their viscosity contrast.

However, reducing the asthenospheric viscosity does help decouple the plate from the underlying mantle. A set of models with uniform lithospheric thickness but varying asthenospheric viscosity confirm this effect (Supplementary Fig. 5), where the mantle drag decreases by a factor of 2.5 (from 5×10^12 Nm^-1 to 2×10^12 Nm^-1) as asthenospheric viscosity drops by a factor of 20 (from 1×10^20 Pa s to 5×10^18 Pa s), consistently below the level of resistance force due to Tibetan GPE. On the other hand, further increasing the asthenospheric viscosity will lead to a saturation in the magnitude of mantle drag (see the relative change in asthenospheric viscosity and mantle drag in Supplementary Fig. 5). This result further confirm the previous finding that, for a given density distribution, buoyancy has a much stronger control on sub-lithospheric stress than mantle viscosity. Consequently, lithospheric thickness variation must contribute additionally to balance the resistance from the Tibetan Plateau.

Previous studies adopting a uniform lithospheric thickness were able to reproduce the present Indian motion, thus raising a question on the respective roles of friction drag, pressure drag, and Tibetan resistance. To explore these forces and the resulting plate motion predictions, we successively reduce the thickness of the continental root. This reduction leads to progressively smaller mantle drags (Supplementary Fig. 6a–c) that eventually fall below the threshold of Tibetan resistance (3×10^12 Nm^-1). With the maximum root depth being 200 km (250 km in the reference case), the Indian plate motion decreases slightly due to a smaller drag of 8.4×10^12 Nm^-1 (Supplementary Fig. 6a). However, truncating the root below 175 km leads to an increase in plate velocity although the drag is even smaller (6.8×10^12 Nm^-1; Supplementary Fig. 6b). This increase of velocity may be attributed to the greater reduction in plate boundary resistance than in mantle drag. If the Indian plate is of a uniform thickness of 100 km (as in a previous study), while the predicted plate motion remains acceptable, the mantle drag falls below the lowest estimate of Tibetan resistance (2.1×10^12 Nm^-1 compared to 3×10^12 Nm^-1; Supplementary Fig. 6c). These results further suggest that the continental root affects the transmission of mantle drag into the overlying lithospheres in two ways: (1) to couple the Indian plate with the underlying asthenosphere, and (2) to maintain a large contact area between Indian and Eurasian plates. Thus, the continental root strongly modulates the balance among the mantle drag, the Indian plate motion, and the northern boundary resistance from the Tibetan Plateau. A recent study shows that the craton is compressed by mantle convection, which is similar to our study except that one side of the cratonic India is sustained by continental collision. Since the mantle drags concentrate on the rooted edge that faces the fast mantle flow, the exact geometry of the root will unlikely affect the overall coupling as long as the flow-normal area of the edge remains similar.

While some researchers suggest that the upwelling plume has a strong impact on mantle drag in this region, our model reveals that the plume arising from the edge of the African LLSVP (Fig. 3f) does not significantly affect dynamic pressure beneath the Indian plate (Fig. 3b). To further examine the effect of mantle upwelling, we perform tests by removing all hot mantle anomalies in the model. Surprisingly, the reproduced Indian plate motion remains similar, and so does the amount of mantle drag (Supplementary Fig. 7). In addition, a recent seismic tomography study shows that the plume is mainly located beneath the Somalian plate, supporting its limited contribution to driving the Indian plate. Therefore, we suggest that sinking slabs play a more crucial role than rising plumes in driving the Indian plate motion.

Following the same model setup that generates the preferred present-day forces (Fig. 4a, b) and kinematics (Fig. 2), we extend the analysis to the geologic past. The calculated mantle drag is consistently larger than the present-day value (~10^13 Nm^-1), with that at 50 Ma being twice as large (Fig. 4c). The calculated drag history closely follows the temporal trends of both predicted and observed plate motions since ~45 Ma (Fig. 4c) when India–Asia collision started in the model (Fig. 1), further confirming mantle drag as the dominant driving force.

The modeled collisional resistance ranges from 4.6×10^12 to 9.7×10^12 Nm^-1, consistent with both present-day GPE estimates (3×10^12 to 8×10^12 Nm^-1) and the magnitude of mantle drag forces. The key finding is the overall consistency of the resistance force magnitude throughout the post-collisional period despite the fluctuation of resistance force at 10 Ma, which is induced by an active boundary-wise lithospheric delamination (Supplementary Fig. 8) temporarily guiding the collisional force downward. This consistency validates our results and agrees with geological evidence for the early Cenozoic uplift of the Tibetan Plateau. The long-term compatibility between modeled forces and independent GPE estimates reinforces the robustness of our findings in capturing the first-order dynamics of the India–Asia collision.

Our estimated plate motion cannot reproduce the much faster Indian motion from 60 to 50 Ma. Both the small magnitude of mantle drag and the assumed single subducting slab (Fig. 3d) during this period implies the operation of other potential driving forces such as those discussed in the following section. The resistance force drops below zero in the pre-collisional snapshots. Because our quantification method does not fully account for the force transmitted through the strong, coherent slab attached to the subducting plate, the negative resistance probably represents a net force contributing to trench advance prior to the continental collision.

Implications on forces driving plate motion and continental collision
The quantification of the Cenozoic mantle drag permits a reevaluation of other potential plate driving forces. Although our model assimilates lateral density variations of oceanic lithosphere following seafloor age, it likely does not fully capture the effect of ridge push due to the flat model surface across the spreading center. However, given its small magnitude of ~3×10^12 Nm^-1 that cannot even overcome the minimum estimate in Tibetan GPE force, ridge push is unlikely to be a major driving force for the Cenozoic Indian motion. We suggest this should be a secondary force that might contribute to explaining the slightly slower predicted plate motion compared to observation (Fig. 4c). Additional results reveal that slab pull from nearby subduction, specifically the Sumatra-Java slab pull (Supplementary Fig. 9) acting through the Indian-Capricorn plate (IND-CAP) boundary (Fig. 1a), should also be minor compared to mantle drag. In practice, the Indian mantle drag and the IND-CAP boundary drag are complementary forces in driving Indian plate motion. With weakened Indian plate boundaries (to minimize surrounding forces), a uniformly thick Indian lithosphere, and a weak asthenosphere (to minimize mantle drag), the Indian plate experiences nearly no motion (Supplementary Fig. 6d). Strengthening the IND-CAP boundary based on this model (Supplementary Fig. 6f) to allow the transmission of the Sumatra slab pull to India only produces a small fraction of the observed plate motion, similar to the scenario where the Indian continental root is thickened by 50 km to increase the lithosphere-asthenosphere coupling (Supplementary Fig. 6e). Since these two models predict similar plate motions, the extra 1×10^12 Nm^-1 mantle drag in the latter model could approximate the amount of IND-CAP boundary drag in the former model, thus confirming its minor role in sustaining the Cenozoic Indian-Tibetan collision.

Our model does not reproduce the faster plate motion (>10 cm yr^-1) observed before ~50 Ma. This may imply the existence of other driving forces besides mantle drag. Relevant hypotheses include a rising plume head, additional slab pull due to double subduction and thinning of Indian lithosphere. Our model does not produce a major plume around this time (Fig. 3d–f), thus not capturing this additional force. The plate reconstruction we adopted does not consider a double subduction scenario either – it instead incorporates a back-arc setting with a southward trench jump. Model tests also show the reduction of lithospheric thickness can increase plate speed with small driving force (Supplementary Fig. 6b). Therefore, we propose that these alternative hypotheses were likely at play during the earliest Cenozoic.

The history of the India–Asia collision provides additional support for the underlying mantle drag as the dominant driving force. The estimated present-day southward resisting force from the Tibetan Plateau (3×10^12 Nm^-1 to 8×10^12 Nm^-1) matches our estimated resistant force from the northern boundary since 40 Ma (Fig. 4c). It also provides a minimum estimate of the required amount of northward Indian plate driving force, substantiating the necessity of mantle drag acted on a thick continental root, since the amplitude of the mantle drag is notably larger than that of other relevant forces as discussed above. The India–Asia collision history was characterized by persistent northward advancement of the Neo-Tethyan subduction zone, which likely started in the Cretaceous when India was initially rifted off from Gondwana. This phenomenon is unique among global plate tectonics where most oceanic trenches tend to roll back. Our results show that a strong northward mantle drag existed even before the onset of collision (e.g., 50 Ma, Fig. 1c). This drag not only maintains the Indian motion but also carries the trench northward, which is a dynamic feature that has not been addressed in most previous studies.

Methods
Governing equations and data-assimilation
The global spherical geodynamic model is based on the code CitcomS. The mantle is divided into 12 caps, each with a high-resolution mesh consisting of 257×257x113 nodes, equivalent to a horizontal (vertical) resolution of 25 km (12 km) near the surface. This resolution allows a proper capture of the lithospheric structure, especially its interaction with the convective mantle. We assume the mantle is incompressible and satisfies the Boussinesq approximation. The model simulates subduction history and mantle evolution through the last 200 Myr. The conservation equations for mass, momentum, energy, and chemical particles are:

∇·u=0

-∇P+∇·(η(∇u+∇Tu))+ρmαΔT+Δρcg=0

∂T/∂t+u·∇T=κΔT

∂C/∂t+u·∇C=0

where u is velocity, P is dynamic pressure, η is viscosity, ρm is the density of the ambient mantle, Δρc is the compositional density anomaly, g is the gravitational acceleration, T is temperature, κ is the thermal conductivity, C is chemical composition. More model-wise parameters are available in the supplementary material (Supplementary Table 1–3).

To reproduce natural, tabular, and asymmetrically subducting slabs that move freely and satisfy observed plate kinematics, we utilize a data-assimilation technique that incorporates a recent plate reconstruction model. The plate motion data from the reconstruction is imposed as surface boundary conditions in the time-dependent model and the bottom boundary is free to slip. The temperature profile of the oceanic lithosphere follows a modified error function following the seafloor age from the reconstruction so that the slab buoyancy is the same as that in the actual Earth, while the continental lithosphere initiates with a steady-state temperature profile (Supplementary Fig. 1a). The core-mantle-boundary temperature is fixed to be 500 °C higher than the ambient mantle. The viscosity of the model depends on depth, temperature, and composition.

The models adopt a layered and smoothed background viscosity profile within 4 depth ranges: 0–44 km, 44–410 km, 410–660 km, and 660–2867 km. Their background viscosities are 10^20 Pa s, 10^20 Pa s, 10^21 Pa s and 3×10^22 Pa s, respectively, and the reference viscosity is set to be 10^21 Pa s (Supplementary Fig. 1b). The temperature dependency of viscosity follows the equation:

η=η0·C·exp(Eη/(T+Tη)-Eη/(Tm+Tη))

where η0 is the background viscosity, C is the compositional viscosity variation, Eη is the activation energy, T, Tη and Tm are the temperature, activation temperature, and ambient mantle temperature, respectively. The variation of viscosity is 10^19–10^23 Pa s in the time-dependent model and 10^18–10^24 Pa s in the instantaneous free-slip models that attempt to reproduce plate motion.

We consider a total of 13 different chemical compositions in the model (Supplementary Fig. 1d), represented by ~1.8 billion tracers. The oceanic lithosphere is composed of a 7-km-thick surface layer that does not have buoyancy anomaly and only mimics the weak and lubricating plate interface near the trench during subduction, a 21-km-thick crustal layer whose buoyancy equals that of a 7-km-thick basaltic crust, and the underlying lithospheric mantle. Both the oceanic basaltic crust and continental lower crust undergoes composition and density changes following the basalt-to-eclogite phase transformation at depths of ~100 km. The continental lithosphere consists of 2 layers of crust and 3 layers of mantle lithosphere. The crust has an average density of ~2800 kg m^-3, with the lower crust being weaker than the upper crust to mitigate the effect of deformation from the imposed plate motion to deeper depths. The continental mantle lithosphere consists of a buoyant upper layer, a neutrally buoyant middle layer, and a dense lower layer relative to the ambient mantle that follows recent inferences. More compositional parameters are available in Supplementary Table 3. The entire continental lithosphere including the crust and the mantle lithosphere within a plate is mechanically coherent because of viscosity cut-off, while the weak zone fully decouples the plate boundary between two plates (Fig. 3c). This allows the stress accumulation in our model to occur primarily within the lithosphere as a whole rather than the crust alone. Thus, our lithospheric resolution is sufficient to capture the first-order force balance. Several recent studies utilizing similar model resolution successfully capture 3D lithospheric stress state, deformation and topography.

Generating plate motion
To reproduce plate motion and stress dynamically, we rerun the model using checkpoints stored at every 10 Myr. These reruns have the same mantle structures as in the sequential model but with a free-slip surface boundary condition. This way, internal mantle dynamics from slabs and plumes control the convective motion and drive deformation of the lithosphere, so that we can estimate the driving force of a given plate. Rigid rotation is removed during each rerun. We further remove the lithospheric net rotation from calculated plate motion for better comparison with observed plate motion in the NNR (no-net-rotation) frame. The mantle flow and stresses from the model with imposed plate motion are largely unchanged in the free-slip reruns, validating the dynamic compatibility of the internal convection with observed surface kinematics. To achieve realistic plateness, we increase the viscosity of the lithosphere to 10^24 Pa s so that each plate moves as a coherent unit with little internal deformation. The viscosity of the oceanic asthenosphere is 10^18 Pa s and plate boundary viscosity is around 10^20 Pa s in order to mimic the effect due to strain-rate weakening. These adjusted shallow viscosity features are not effectively constrained in the sequential simulation when imposed surface kinematics override the dynamics of these features. However, a proper representation of them in the free-slip reruns is crucial given their high sensitivity to plate motion and mantle dynamics.

Quantification of mantle drag and boundary resistance
Mantle drag is a combination of friction drag and pressure drag. Friction drag is essentially the shear stress applied near the bottom of the lithosphere due to asthenospheric flow. The friction drag along plate motion direction follow the form:

Ff=∬(τrθupθ+τrϕupϕ)dA

where Ff is friction drag, up is plate motion velocity, τ is the stress, r, θ, and ϕ are the radial, colatitudinal, and longitudinal direction, respectively, and A is the horizontal area. This value is positive for driving the plate motion and negative for resisting the plate motion. The stress used here is the maximum absolute value over depths near a viscosity-defined LAB (viscosity 10 times that of the background viscosity), from 0.75 times of this depth to 410 km.

Pressure drag is due to the blocking of horizontal mantle flow by the lithospheric step. This reflects the pressure gradient acting on the lithospheric step. Thus, computing pressure drag requires a 3D integration of the pressure gradient over the lithospheric step direction as:

Fp=∭ωup·nup(∇P·n)dV

where Fp is pressure drag, ω is a weighting function from the lithospheric step, n is the lithospheric step normal vector (pointing to the lithospheric interior) and V is the volume of interest. Again, the positive value represents a driving force and vice versa. Before computing the integration, the algorithm goes through each depth of interest to find the distribution of the lithosphere (set as 1) and asthenosphere (set as 0) horizontally and acquire its gradient map (the lithospheric step), n is the gradient unit vector and ω is the scalar field of gradient amplitude normalized by the maximum value in the region of interest, ranging from 0 to 1.

During the calculation of plate driving force, we neglect the area within 300 km from the plate boundary to avoid potential boundary effect. To calculate the total force acting on the plate, we first convert the force, including friction drag and pressure drag, from each unit into torque, sum the torques together, and then convert the total torque back into an equivalent force using a reference vector. This reference vector points towards the center of the plate and has a length equal to the Earth’s radius (6371 km). For comparison to traditional estimation of slab pull and ridge push, all the mantle drag forces (in the unit of N) are divided by the square root of plate area in each time to give the unit of Nm^-1.

Plate boundary resistant force is computed similar to pressure drag. The northern resistant force reveals as a dynamic pressure increase towards the Indo–Eurasia collision zone. Thus, we track the pressure force within 300 km and above a given latitude from the plate boundary to quantify northern boundary resistant force. The resistant force has different sign compared to driving force when projected to plate motion.

Anisotropy calculation
We simulate the flow-induced seismic anisotropy based on the FORTRAN code DrexS that considers the effect of dislocation creep, dynamic recrystallization, and grain boundary migration, similar to previous works. The simulation initializes with randomly oriented mantle mineral aggregate that assumes 70% A-type olivine and 30% enstatite. These mineral aggregates advect and realign following the mantle flow from 20 Ma to the present day based on the free-slip models. This produces the full elastic tensors of the aggregates and, since the upper mantle is dominated by transverse isotropy, we use the symmetry axis of the transverse isotropy (TI axis) to represent the lattice-preferred orientation.

Test on free-slip and free-surface model
A regional 2D thermochemical model simulates continental collision driven by detached, sinking slabs in the lower mantle. We use the open-source software ASPECT to solve the governing equations of conservation of mass, momentum, energy, and composition. An unstructured mesh resolves the 6000×2900 km numerical model domain with a grid resolution of 2.9×2.8 km maximum per element. Both experiments use visco-plastic non-Newtonian rheology laws. Both models run for a 0.4-Myr duration. The model includes creeping flow with thermal and chemical buoyancy forces, accounts for eclogite, olivine-spinel, and spinel-perovskite phase transitions, and considers the effects of adiabatic, shear, and radioactive heating.

Two numerical experiments were conducted with free-slip and free-surface boundary conditions for the upper surface, while all other conditions remain identical between the two models. The model prescribes two 2400-km-wide continental plates, separated by a 12-km-wide weak zone in the middle, in the central top of the model. The gap between the plate edge and the model boundary is 600 km. Two vertical weak zones are set in the continental plates to allow the plate to move and to be subducted. The mantle flow that induces active continental collision. The slab has an effective density of 3500 kg m^-3.

The model continental crust has a thickness of 41 km and consists of three layers. Felsic rheology of wet quartzite represent the upper (10 km thick) and middle (15 km thick) crust, whereas mafic rheology of plagioclase represents the lower crust (16 km thick). Anhydrous peridotite with the rheology of dry olivine represent the mantle lithosphere.

The half-space cooling model with an age of 200 Myr defines the initial thermal structure of the continental lithosphere within 1500 km to the continental collision zone. To effectively simulate the effect of lithospheric step, we set the initial thermal structure for the further 900 km of the continental lithosphere to be equivalent to an age of 80 Myr. The initial temperature in the background mantle is 1300 °C plus an adiabatic gradient of 0.3 °C km^-1.

In both models, all mechanical boundary conditions are free to slip, except for the model with a free top surface. The Arbitrary Lagrangian–Eulerian (ALE) method tracks the top surface for the model with a free-surface boundary condition. Erosion and sedimentation, according to the hillslope transport diffusion method with a coefficient of 10^-8 m^2 s^-1, cause the evolution of this upper boundary.

The density of rocks varies with temperature T according to the equation:

ρ(T)=ρ0(1-α(T-T0))

where ρ0 is the reference density at T0 = 298 K, and α is the thermal expansion coefficient.

The model considers the transformation of olivine to wadsleyite and ringwoodite and to bridgmanite in the mantle. An increase in the density of the crust to 3550 kg m^-3 considers the eclogitization of the mafic lower continental crust. The model does not consider melting or hydration processes. The ρ0 for upper crust, middle crust, lower crust, and mantle are 2700, 2800, 3000 and 3300 kg m^-3, respectively.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "I enjoyed reading this manuscript. It is well written and well-illustrated. The videos provided in supplement are quite helpful. The study addresses an important question in geodynamics. The geodynamical mode is very advanced and described in a accessible way with the right level of details. The authors suggest that the plate motion of India is driven primarily by drag imparted by a large scale mantle convection cell. While the basic idea is similar to the ‘conveyor belt’ idea of Becker and Faccenna (2011) the origin of the cell is very different and the authors make a good case that such a large scale convection cell might indeed exist. So I think the manuscript can be considered for publication. I do however have a number of concerns. I am returning an annotated file with detailed comments and suggestions. Most of them are editorial in nature. I only report my main concerns here.\n- My main concern with study regards a possible bias introduced by the fact that CitcomS doesn’t allow to account for the free surface and variations of GPE. The author find a resistive stress along the Himalaya that is of the order of the resistive stress derived from the crustal thickness of Tibet and report this finding as a sanity check. I am unconvinced that this is a good validation that the GPE is somehow rendered in a satisfying way. I note that the purple curve in Figure 4d implies a resistive force twice the present value at 10Ma. I am sure that the authors would be hesitant to give them any geological significance to such a variation. At face value would mean a Tibet half as high as at present by then, while most geological observations suggest that Tibet rose to an elevation close to present quite early in the collision. It therefore probable that the variations of that curve rather show the uncertainties in the estimate of this resisting force. In that case the fact that the estimated value at present matches approximately the value estimated from the crustal thickness of Tibet seems coincidental to me.\n- The authors don’t plot the value of the resisting force before 30Ma because of slab break-off happening at that time. This is exactly the kind of information that I would like to see reported and discussed. In principle slab-break off must have resulted in a sudden drop of slab pull. Because plate motion is assimilated in the model, there is not slow down of India happening in the model. I guess that there must then be a significant drop of the resisting force along the Himalaya front. If so would that make any sense? I’d like the authors to report and discuss what is happening then\n- It makes little sense to quantify the effect of mantle drag in terms of a force per unit length. Reporting a force per unit length makes sense only for plate boundary forces. The authors should report the force and the torque obtained by integrating the sub-lithospheric shear stress over the plate. They should also discuss the force (and torque?) budget taking all sources of stress applied to the India plate. The stresses transmitted across the boundaries with the Arabia and the Capricon plates should be taken into account in that budget.\n- The plate boundary forces (slab pull and ridge push) reported in the study are comparable to those we reported in Copley et al (2010) (when expressed in force per unit length), and so is the resistive force along the Himalayan front. We however concluded that basal drag had to be resistive. I wonder why with similar forces per unit length this study requires such a large mantle drag. I don’t contest the finding but I don’t clearly see the explanation.\n- I am puzzled by the apparent inconsistency between the rather steady plate velocity predicted by the model in Figure 4 and the nearly null velocity seen at 40Ma in Figure 3. The authors need to specify the reference frame used to estimate velocities (NNR, HS, other?) in these figures.\n- The fit between the predicted and observed Transverse Isotropy is a nice result but note that it says that plate motion is consistent in direction but the sense of shear could be either resisting or driving in principle.\n- As I recall Bird 1998 also concludes, based on a global kinematic analysis ,that plates are driven by mantle drag and make a comparison with plate motion. It would probably be appropriate to check and discuss this analysis.", "reviewer_2": "This is an interesting study in which the authors use a classical geodynamic mantle convection code, CitcomS, to make the case that large scale mantle flow is what keeps the Tibetan Plateau up.\nThe idea that a large scale mantle wind drives plates is not new; Faccenna and Becker (2010 Nature) have popularized this for the Mediterranean and applied it later to the India-Asia collision zone in Becker and Faccenna 2011 (using the same code). The results in the current study largely confirm this and also demonstrate that slabs in the lower mantle play a key role in driving the plates. The study is an improvement over previous work because 1) a novel data assimilation method is employed that takes plate reconstructions into account; 2) a more serious effort is done to quantify the forces that occur in the simulations.\nThis is clearly of interest to a broad audience, and the codes and data assimilation method have been used in many previous studies.\nYet, I do have two major concerns, that can be addressed with additional simulations, as well as a few more minor remarks:\n1) A key aspect of this paper is that the gravitational potential energy of the Tibetan plateau, which is caused by the difference in elevation of the plateau and the surrounding lower areas, is very significant and on the order of 3-9x1012 N/m as discussed in the main text. The goal of the paper is to quantify mechanisms that are able to provide sufficient forces and the authors do a good job in that. Yet, a key limitation of their models is that it does not include a free surface and as a result, the gravitational potential energy (the resistive force) is not taken into account. The authors are aware of this limitation and point out that a free surface is not being used (lines 160-163). Yet, they only argue on the basis of the forces that are estimated from natural observations, while using the driving forces from their model without topography. In models with an actual free surface, the magnitude of these driving forces may actually be quite different from those of free slip models; I am not aware of a study that has quantified this error and therefore I cannot judge whether the key conclusions of this study will hold. If you implement a setup with a changing crustal thickness and a free slip top boundary, the crust wants to move upwards due to a density instability that causes additional forces. As such\nWhereas I realize that it is unfair to ask the authors to change CitcomS to take a real free surface into account, one could estimate the error by performing 2D or 3D idealized simulations with other open-source geodynamic codes (such as ASPECT, pTatin3D, or LaMEM) that do have a self-consistent free surface implemented. It would be rather straightforward to do a similar force balance analysis for a setup with a craton that has roots of varying thickness and a large scale mantle flow, and uses either free slip or free surface upper boundary conditions. D oing\n2) The resolution of the models is rather coarse with only 113 nodes in the vertical direction, and a vertical resolution of around 12 km in the crust. This is hardly 4 elements over the whole crust which is on the lower bound of what is required to even resolve this (given the linear elements that CitcomS uses). It is clearly crucial that the authors perform resolution sensitivity tests to demonstrate that their results are robust.\nMinor remarks:\n- For how long have the forces been sufficient large to keep the Tibetan plateau high?\n- Please use SI units throughout (2800 kg/m3 instead of 2.8 g/cm3)\n- Please convert viscosity values on Fig. 3c into dimensional units (Pas)\n- Please add a table with all material parameters employed. The statement “ more model information is available in recent studies using similar methods” is very unsatisfying as it doesn’t allow me to reproduce what was actually done."}, "author_response": "We thank the two reviewers for their constructive comments. We have now systematically revised the manuscripts based on these reviews. The original comments are in black font and our responses in blue. In response, we have made substantial revisions to the manuscript. The main changes are summarized as follows:\nWe have included a comprehensive description justifying the use of the free-slip model in understanding plate force balance, considering the influence of gravitational potential energy arising from topography (new Supplementary Fig. 10). Additional models were run to compare the force profile across the collisional boundary, demonstrating that the force distribution between a free-slip and a free-surface model is highly similar.\nThe complete resistance force history (prior to 30 Ma) is now presented to illustrate the change occurring during continental collision (Fig. 4c). Further analysis of the snapshot at ~10 Ma revealed that the anomalous drop in resistance force is attributed to a lithospheric delamination event on the Eurasian side, which diverts the resistance force downward. This finding supports the overall consistency of the resistance force.\nWe have refined the language of the manuscript, enhanced the introduction, ensured consistent use of units, and improved the reasoning logic throughout. We have also provided a throughout elaboration on sensitivity concern raised by the second reviewer. Additionally, we have included three tables (Supplementary Tables 1, 2, and 3) to provide complete model parameters. Furthermore, we have expanded the discussion to highlight the advantages and potential areas for further improvement in our current geodynamic models.\nWe hope the additional work addresses all the reviewers’ concerns. The following are detailed one-to-one responses.\nReviewer #1\nI enjoyed reading this manuscript. It is well written and well-illustrated. The videos provided in supplement are quite helpful. The study addresses an important question in geodynamics. The geodynamical mode is very advanced and described in a accessible way with the right level of details. The authors suggest that the plate motion of India is driven primarily by drag imparted by a large scale mantle convection cell. While the basic idea is similar to the ‘conveyor belt’ idea of Becker and Faccenna (2011) the origin of the cell is very different and the authors make a good case that such a large scale convection cell might indeed exist. So I think the manuscript can be considered for publication.\nAuthor response: Thank you very much for the encouraging comment.\nI do however have a number of concerns. I am returning an annotated file with detailed comments and suggestions. Most of them are editorial in nature. I only report my main concerns here.\nAuthor response: Thanks for the detailed suggestions, and we have responded to all the annotations after addressing the main concerns.\n- My main concern with study regards a possible bias introduced by the fact that CitcomS doesn’t allow to account for the free surface and variations of GPE. The author find a resistive stress along the Himalaya that is of the order of the resistive stress derived from the crustal thickness of Tibet and report this finding as a sanity check. I am unconvinced that this is a good validation that the GPE is somehow rendered in a satisfying way.\nAuthor response: Thank you for raising this important point about the potential limitations of using a free-slip surface boundary condition in CitcomS. We agree that a free-slip model has some difference from a free surface model in capturing the gravitational potential energy (GPE). However, we believe that our key findings regarding the plate-wise force balance, particularly the magnitude of mantle drag and the resistance related to GPE, remain robust for the following reasons:\n1. Previous studies have shown that the difference in mantle flow between free-slip and free-surface models is relatively minor (Zhong et al., 1996). This suggests that the resulting plate-driving forces, especially due to mantle drag, should be similar between the two types of models. Mechanically, a free-slip boundary, although with a flat surface, stores a finite vertical stress component that is equivalent to the stress as the same depth of a free-slip model; this means their respective GPE distribution is also similar. Indeed, codes like CitcomS compute surface topography by converting this sub-surface stress into a topography that very much resemble the real case, as seen from the figure below (Review Fig. 1) that we just published.\nReview Figure 1. Predicted topographic evolution of North America based on a data-assimilation model from CitcomS (Li et al., ESR, 2023). The modeled topography resembles the real Earth topography both in the amplitude (>5 km) and geological observations (red arrows at 120 Ma mark paleo current direction, and red lines at 80 Ma represent paleo shorelines).\n2. Some earlier studies show that free-slip models can generate more short-wavelength oscillations in topography compared to free-surface models (e.g., Fig. 4 in Crameri et al., 2017), but the longer-wavelength topographic features are generally consistent. Zhong et al. (1996) show that the smoothing in topography in free-surface case is due to the temporal retard from surface relaxation. In the Crameri et al. (2017) study, the slab geometry is somewhat artificially defined, i.e., not naturally produced by the model. We suspect, if the slab was generated more naturally (thus better trained by the boundary condition), the calculated (from sub-surface stress) surface topography would be more consistent with the free-surface case, as shown below in our new test. In fact, even assuming such short-wavelength topographic differences do exist, they would have a small effect on the overall plate force distribution. Since our analysis focuses on the plate-scale force balance rather than local topographic variations, these short-wavelength differences are unlikely to significantly impact our overall conclusions. On top of this, such short-wavelength oscillation is also related to the rough surface viscosity variation in similar spatial scale (e.g., Fig. 4a,b,d,e in Crameri et al., 2017, likely due to plastic deformation that amplifies the intrinsic rheological inconsistency with the arbitrarily defined slab geometry), which is not the case in our model (Fig. 3c).\n3. To further test the robustness of our results, we conducted an idealized numerical experiment using the geodynamic modeling software ASPECT. We ran both free-slip and free-surface models of mantle-drag-driven continental collision over 0.4 Myr. As expected, the free-surface model produced a smoother topographic profile while the free-slip model exhibited more short-wavelength features (Review Fig. 2b). Crucially, however, the vertically integrated lithospheric pressure, which represents the accumulation of driving and resisting forces, was very similar between the two models (Review Fig. 2a). Additionally, the mantle flow patterns and plate velocities were nearly identical (Review Fig. 2c, d). All these confirm our previous inference and support the force analysis in a free-slip model.\nThese results give us confidence that, although our CitcomS models do not explicitly include a free surface, they still capture the first-order effects of GPE on the lithospheric stress state and plate-driving forces. The strong similarity in vertically integrated pressure between the free-slip and free-surface ASPECT models supports the validity of the resistive force we infer along the Himalaya plate boundary. We have added in line 319-336 “Previous research suggests that mantle models based on a free-slip… as it captures the essential features of the stress distribution” to elaborate on this in the main text, and included Review Fig. 2 as our new Supplementary Fig. 10.\nWe acknowledge that future work is needed to fully establish the force variation between a free-slip and free-surface model being minor. In this sense, the additional ASPECT models are beyond the scope of the current manuscript and will be developed into a separate publication.\nReview Figure 2. parison of results between the free slip and free surface models. (a) Depth integrated lithospheric dynamic pressure and (b) topography for both the free-slip and free-surface models. Lithospheric dynamic pressure is integrated vertically from the surface to 1000°C isotherm within the lithosphere, which is roughly equivalent to the 1e21 viscosity boundary. It directly reflects the lithospheric force accumulation along the horizontal direction. (c) Free-slip and (d) free-surface models showing viscosity and velocity.\nI note that the purple curve in Figure 4d implies a resistive force twice the present value at 10Ma. I am sure that the authors would be hesitant to give them any geological significance to such a variation. At face value would mean a Tibet half as high as at present by then, while most geological observations suggest that Tibet rose to an elevation close to present quite early in the collision. It therefore probable that the variations of that curve rather show the uncertainties in the estimate of this resisting force. In that case the fact that the estimated value at present matches approximately the value estimated from the crustal thickness of Tibet seems coincidental to me.\nAuthor response: Thank you for pointing out the potential implication of the resistance force fluctuation. We have conducted further analysis around the 10 Ma snapshot and noticed a large-scale Eurasian lithospheric delamination event occurring along the collision boundary. This delamination is likely the main cause of the reduced resistance, as the force is diverted downward and acts less horizontally into the collision zone. We have added a new supplementary fig. 8 to demonstrate this process and included the following statement in line 230-233: \"The key finding is … despite the fluctuation of resistance force in 10Ma, which is induced by an active boundary-wise lithospheric delamination (Supplementary Fig. 8) guiding collisional force downward. \"\nHowever, since our model does not allow free-deformation within the continental interior (due to imposed rigid-plate-style motion), this delamination event may not reproduce precisely a realistic event, although delamination is expected to occur at some point during the continental collision process. We prefer to relay this topic to future research, when further implementing freely deformable plate interiors. Therefore, we consider the temporal fluctuation in resistance force caused by delamination to be plausible, even if its timing may not be fully realistic.\nOn top of this, we emphasize that despite these short-term fluctuations, the key point is that the modeled resistance force remains within the range of 3×1012 to 8×1012 N/m (derived from GPE estimates) throughout Cenozoic (see updated Fig. 4). This indicates that our modeled forces are consistently of sufficient magnitude to support an elevated Tibetan Plateau over the entire post-collisional period (~40 Ma to present), in agreement with geological evidence for an early Tibetan uplift after collision.\nWe argue that the fact that these values remain compatible over tens of millions of years is a robust result and not merely coincidental. This long-term consistency, rather than the agreement at single timestep, is the key evidence that our models capture the salient dynamics of the India-Asia collision.\nWe have also updated the manuscript to clarify these points and avoid placing undue emphasis on short-term variations in line 228-236 “The modeled collisional resistance … first-order dynamics of the India-Asia collision.” We hope this addresses the reviewer's concerns.\n- The authors don’t plot the value of the resisting force before 30Ma because of slab break-off happening at that time. This is exactly the kind of information that I would like to see reported and discussed.\nAuthor response: Thank you for raising this point. We have now updated Fig. 4 to include the resisting force values before 30 Ma, providing a more complete picture of the force variation over time.\nHowever, it's important to note that our method for computing boundary resistance only reflects the total force acting at the plate boundary when there is no coherent slab attached to the subducting plate. When a slab is present, the force transmitted through the slab is not fully accounted for in our calculation.\nThe negative resistance force observed before 40 Ma likely reflects the dynamics of trench advance. This negative force may be related to the dynamics inducing the northward motion of the Neo-Tethyan subduction zone prior to the continental collision.\nThese pre-30 Ma force variations provide additional context for the evolution of the India-Asia collision. However, our main conclusions regarding the magnitude and consistency of the collisional resistance during the post-collisional period (40 Ma to present) remain unchanged.\nWe added “Because our quantification method does not … trench advance prior to the continental collision.” in line 241-243 to clarify the implication for the result prior to 40 Ma.\nIn principle slab-break off must have resulted in a sudden drop of slab pull. Because plate motion is assimilated in the model, there is not slow down of India happening in the model. I guess that there must then be a significant drop of the resisting force along the Himalaya front. If so would that make any sense? I’d like the authors to report and discuss what is happening then\nAuthor response: Thank you for raising this intriguing question about the potential effects of slab break-off on the force balance and plate motion. While a detailed analysis of the full force balance, including slab pull, during the slab break-off event is indeed an interesting topic, it is beyond the scope of the current paper.\nIn our model, the Neo-Tethyan slab is in a rollover geometry at the time of break-off. This configuration likely results in a less efficient transfer of the gravitational force from the dense slab to the plate, compared to a slab in a rollback geometry. Consequently, the reduction in slab pull during break-off may be more gradual and less pronounced than in a scenario with a more favorably oriented slab.\nMoreover, our model suggests that the mantle drag driving force initiated before the collision (for example, at 50 Ma in Fig. 4c). The combination of a less efficient slab pull and the early onset of mantle drag likely results in a gradual transition from a slab-pull-driven to a mantle-drag-driven plate motion, rather than a sudden drop.\nHowever, it is important to note that our model does not fully reproduce the rapid Indian plate motion observed prior to the collision. As a result, any discussion of the complete force balance before slab break-off may not provide a realistic representation of the tectonic evolution during this period.\n- It makes little sense to quantify the effect of mantle drag in terms of a force per unit length. Reporting a force per unit length makes sense only for plate boundary forces. The authors should report the force and the torque obtained by integrating the sub-lithospheric shear stress over the plate. They should also discuss the force (and torque?) budget taking all sources of stress applied to the India plate. The stresses transmitted across the boundaries with the Arabia and the Capricon plates should be taken into account in that budget.\nAuthor response: We agree that mantle drag, being a distributed force acting over the base of the lithosphere, is fundamentally different from plate boundary forces like slab pull or ridge push. In our study quantifying the total force, we initially calculate the mantle drag force as a torque by integrating the sub-lithospheric shear stress over the entire plate. This approach provides a comprehensive measure of the total force acting on the plate due to mantle drag. We then convert this torque into an equivalent force per unit length by dividing it by the square root of the plate area. We updated a clarification in line 464-468 “To calculate the total force acting on the plate, we first convert the force, including friction drag and pressure drag, from each unit into torque … This reference vector points towards the center of the plate and has a length equal to the Earth's radius (6371 km)”\nWe chose to report the mantle drag force in terms of force per unit length primarily to facilitate comparison with other well-known plate driving forces, such as slab pull, ridge push, and forces associated with gravitational potential energy (GPE) variations along plate boundaries. These forces are commonly reported in the literature as force per unit length, and adopting the same convention for mantle drag allows for a more straightforward comparison of their relative magnitudes and contributions to the overall force balance.\nHowever, we acknowledge that this conversion may not fully capture the distributed nature of the mantle drag force. To provide a more complete picture of the mantle drag force distribution, we also present the force in terms of equivalent stress in Fig. 4. While friction drag itself is acquired from the maximum shear stress in the vertical axis, the pressure drag is quantified from pressure gradient, which can only be equivalently converted to the unit of a stress (yet not a real stress). This representation allows readers to visualize the spatial variability of the mantle drag force across the plate.\nWe recognize that directly reporting the torque would avoid any potential confusion arising from the conversion to force per unit length. However, we believe that the benefits of enabling direct comparisons with other plate driving forces outweigh the potential drawbacks, especially when supplemented with the stress distribution map.\nThey should also discuss the force (and torque?) budget taking all sources of stress applied to the India plate. The stresses transmitted across the boundaries with the Arabia and the Capricon plates should be taken into account in that budget.\nAuthor response: Thank you for raising this important point about considering the stresses transmitted across the boundaries with the Arabia and Capricorn plates in the force budget. We agree that understanding the full force budget is crucial, and we appreciate the suggestion to include these boundary stresses.\nWe did not report the forces in these plate boundaries for two main reasons:\n1. When quantifying the resistance force in a manner similar to the northern boundary, we found that the forces in the Arabia and Capricorn plate boundaries were significantly smaller than the resistance in the north and the mantle drag. Using the same color scale as the northern resistance, the boundary forces in the Arabian and Capricorn plate boundaries would be mostly invisible.\n2. Our quantification only considers normal pressure forces acting on the plate boundary and projects them onto the plate motion direction. This means that shear stresses may not be fully accounted for. For example, the Capricorn plate motion is nearly parallel to the Indian plate motion, suggesting that shear stress may be acting on the boundary instead of normal pressure force, [Redacted].\nNevertheless, we acknowledge that confirming the forces transmitted through the Capricorn and potentially Arabian plate boundaries are mostly negligible is crucial for our conclusion. This directly leads to our reasoning in line 253-266, demonstrating the minor role of these forces by fully weakening the IND-CAP plate boundary, which potentially eliminates all shear stress transferred through that boundary.\n[Redacted]\n[Figure Redacted]\n- The plate boundary forces (slab pull and ridge push) reported in the study are comparable to those we reported in Copley et al (2010) (when expressed in force per unit length), and so is the resistive force along the Himalayan front. We however concluded that basal drag had to be resistive. I wonder why with similar forces per unit length this study requires such a large mantle drag. I don’t contest the finding but I don’t clearly see the explanation.\nAuthor response: Thank you for raising this point. In Copley et al. (2010), the Indian plate, as defined in our study, along with the Capricorn and Australian plates, were collectively referred to as the \"Indian plate\" in the present-day context. Therefore, the slab pull force reported in Copley et al. (2010) likely originates from the slabs attached to the Capricorn and Australian plates, which are not included in our calculation.\nWe agree that if we consider the Australian and Capricorn plates together with the Indian plate as a single \"Indian plate,\" the dominant force may indeed be slab pull, since the fast mantle flow only appears beneath the Indian plate, but the Australian and Capricorn plates combined take up more area than the Indian plate, thus resulting in a smaller average basal drag. However, our test of the Indian-Capricorn diffusive plate boundary, as mentioned in line 253-266, demonstrates that totally weakening this boundary does not change the Indian plate motion. This suggests that although the Indian plate motion is similar to the Australian and Capricorn plate motions, it is mostly self-sustained by its own mantle drag driving force.\nIn other words, according to our research, the Indian post-collisional motion is not sustained by the slab pull force transmitted from the Capricorn and Australian plates. The apparent similarity in the plate motions of the Capricorn, Australian, and Indian plates may be coincidental, and treating them as a single plate could obscure the more complex dynamic mechanisms at play. [Redacted]\n- I am puzzled by the apparent inconsistency between the rather steady plate velocity predicted by the model in Figure 4 and the nearly null velocity seen at 40Ma in Figure 3. The authors need to specify the reference frame used to estimate velocities (NNR, HS, other?) in these figures.\nAuthor response: Thank you for pointing out the apparent inconsistency between the plate velocities shown in Figures 3 and 4. We apologize for any confusion this may have caused and appreciate you bringing this to our attention. The discrepancy arises from the use of different reference frames in these figures. In Figure 2, we used the No-Net-Rotation (NNR) reference frame, while in Figure 3, we used a mixed reference frame based on Müller et al. (2016), which is a mixed reference frame that includes hot spot information. When combined with net rotation, the plate motion direction becomes nearly eastward, resulting in the very slow velocity shown in a latitudinal cross-section. To clarify this, we have added the following sentences to the Figure 3 caption:\n\"The velocities in b and c are from free-slip models and follow a no-rigid-rotation reference frame.\"\n\"The velocities here are from the imposed plate motion model and follow the reference frame used in the plate reconstruction2.\"\nThank you again for your careful review and for bringing this inconsistency to our attention.\n- The fit between the predicted and observed Transverse Isotropy is a nice result but note that it says that plate motion is consistent in direction but the sense of shear could be either resisting or driving in principle.\nAuthor response: We agree with the reviewer that while the fit between the predicted and observed Transverse Isotropy is consistent, it does not necessarily indicate whether the shear is representing driving or resisting force to the plate motion. The fit is a necessary but not sufficient condition for mantle flow to be the driving force, providing a partial validation of the model outcome. We have changed the sentence in line 89-90 into “… seismic anisotropy (Fig. 2d) further supports the directionality of this fast mantle flow” to relax the argument.\n- As I recall Bird 1998 also concludes, based on a global kinematic analysis, that plates are driven by mantle drag and make a comparison with plate motion. It would probably be appropriate to check and discuss this analysis.\nAuthor response: Thank you for highlighting the work by Bird (1998), which concludes that the Indian plate is partially driven by mantle drag based on a global kinematic analysis, and that the mantle drag likely acts primarily on the continental area. To provide a bit more information, Bird (1998) approximates mantle flow velocity by using a scaled-up plate motion velocity, which, while kinematically showing that mantle drag is likely contributing to the driving force, does not reveal the true dynamics of such mantle flow.\nTherefore, it is encouraging to see that our results confirm Bird's (1998) proposal and take a step further by discussing the origin, dynamics, and temporal evolution of the mantle flow beneath the Indian plate. We have added discussions of Bird's (1998) work (new ref. 25) in line 138-147 and 161-163 of our manuscript highlighting the consistency between our findings and their conclusions.\nComments in the pdf file:\n(Line numbers from reviewer’s comment are from the old manuscript while author response uses new line numbers in the updated manuscript.)\nLine 21. This is not a natural unit to quantify mantle drag. Mantle drag is better quantified in term of the basal shear stress. You can also integrate over the plate area to get the total force but there is no natural length to convert mantle drag into a force per unity length. (the mean width of the Indian plate measured perpendicular to the plate motion relative to the upper mantle is not a very natural length...).\nAuthor response: Thanks for bringing this point, we have addressed this concern in the main comment section.\nLine 23. I find that result puzzling. Ideally you would like to enforce the GPE based on the known topography and have a real free surface. I understand that this is not possible in a global geodynamical model. However the result seems coincidental to me.\nAuthor response: Thanks for bringing this point, we have addressed this concern in the main comment section.\nLine 30. It is fine to quote recent plate kinematic models but this fact has been known for a long time. It would be appropriate to quote one of the seminal articles which reported it (possibly Molnar and Tapponnier, 1975).\nAuthor response: Thank you for bringing up this important seminal article. We agree with your point and have added the mentioned seminal work to the introduction.\nLine 33. This is a very vague statement.\nAuthor response: We reorganized this statement to the end of the paragraph in line 46-48 “Thus, the post-collisional Indian motion cannot be adequately explained by the above mentioned plate driving forces8, leading to the intriguing question of what the dominant driving mechanism and its associated origin are.” We hope the scientific question now is better presented.\nI had not seen reference #3. I just looked at it to have a clearer idea of the 'questions' you are referring to. It seems to me it is an opinion piece with no particular question highlighted. It doesn't help define the motivation for your work.\nAuthor response: Thank you for the comment. Following your previous comment, we have reorganized the paragraph. The revised paragraph now more clearly highlights the key knowledge gap and the specific question we aim to address in this study without relying on the previously mentioned reference.\nPlease state more clearly what is the intriguing question you want to address. I have had interest in this topic for a long time and my opinion is that we understand rather well the different factors at play but quantification remains a challenge. For example, some analysis have argued that 'basal drag' is resisting the India plate motion (this is one point we make in Copley et al), other have argued that it is a driving force (e.g., Becker and Faccenna, 2011).\nAuthor response: We agree that to some degree, we know what possible factors can potentially be at play. And it is precisely the quantification that causes the intriguing question to persist: to what extent do these factors contribute, and could we attribute any of them as a dominant driving force? Furthermore, what are the origins of these different factors? For example, if mantle drag is a significant driving force, which anomaly drives the fast mantle flow? We hope the updated statement more clearly conveys this concept.\nLine 35. Might be more appropriate to quote Lyon-caen and Molnar, 1988 here rather than ref 4. This is where the GPE due to the thickened crust in Tibet was first evaluated and more recent estimates are not much different (as I recall the value reported in ref 4 is actually taken from this previous study). Note that the central and provocative idea proposed in ref 4 is that delamination of the lithosphere and uplift of Tibet would have occurred in the mid-Cenozoic in association with strengthening of the Monsoon. This is not a topic that you are discussing and it is somewhat outdated anyway (the mantle lid is still present beneath most of Tibet as shown by surface wave tomography, and there is not evidence for a sudden uplift of Tibet in the mid Cenozoic).\nAuthor response: Thank you for the suggestion. We agree and have added the Molnar and Lyon-Caen, 1988 paper to the citation. This paper lays the foundation for estimating the force per unit length and provides an initial estimate of 6.9x1012 N m-1. However, while it is true that the central and provocative idea proposed in Molnar et al. 1993 (old ref. 4, new ref. 5) is indeed not a topic within our discussion, we have also retained the Molnar et al. 1993 (old ref.4, new ref. 5) paper because it offers a more thorough estimation with a range of 8±2x1012 N m-1, which serves as a useful upper limit that we wish to include.\nLine 36. I don’t get what you mean here. Ref 7 is a very important paper, probably the first attempt at quantifying plate boundary forces based on kinematics. So it is fine to quote but the statement makes no clear sense to me.\nAuthor response: In the highly influential paper by Forsyth and Uyeda, 1975 (old ref. 7, new ref. 8), the authors discuss the various forces contributing to plate force balance, including different types of resisting forces that counteract the driving forces. Collisional resistance is one of the resisting forces considered in their analysis. In this context, a higher driving force is required to overcome not only the collisional resistance but also other types of resisting forces to maintain the observed plate motion. This is the reason we cite this reference here, as it provides a framework for understanding the balance between driving and resisting forces in plate tectonics.\nLine 38. Dismissing the importance of slab pull is too casual.\nAuthor response: The statement in the old line 38 (new line 43) aims to clarify that direct oceanic slab pull along the northern plate boundary of the Indian plate is not a significant driving force in our study, as we have defined the Indian plate to be separated from (via a weak boundary) the Capricorn and Australian plates. However, we acknowledge that slab pull from the Capricorn and Australian plates can potentially be transmitted to the Indian plate and act as a driving force. This possibility is discussed in the following paragraph and is further explored through model tests and analysis in lines 253-266. Our findings suggest that these transmitted slab pull forces play a minor role in driving the India-Eurasia collision.\nLine 41. I disagree with this statement. In our analysis we found that the force per unit length associate with subduction and ridges where very close. If you compare plate rate (relative to HS or NNR) with the fraction of the plate boundaries where a driving force is applied due to either subduction or oceanic rifting the correlation is very compelling (figure 6 in Copley et al). I think many studies are just underestimate the effect of the GPE at mid oceanic ridges. I have followed closely recent work on this issue though.\nAuthor response: As we have addressed a similar comment earlier, we would like to reiterate that we agree with the analysis from Copley et al., 2010, regarding the importance of slab pull as a major plate driving force when considering the Indian plate to include the Capricorn and Australian plates, assuming the slab pull from the latter can readily transmit to the Indian plate. However, in our study, we demonstrate that the plate motion of the Indian plate, as we have defined it, is primarily maintained by mantle drag, while the Capricorn plate and possibly the Australian plate are predominantly driven by its own slab pull. [Redacted] Therefore, our findings suggest that these plates are likely in their respective force balance states, providing a new insight that does not contradict the analysis from Copley et al., 2010.\nLine 44. This is true only if the crust is eclogitized and retain enough strength. It is not clear that substantial subduction of continental crust has occurred. Some geological evidence suggest that slab break-off occurred relatively early (De Sigoyer et al, 2000; Guillot et al, 2003). The formation and exhumation of eclogites occurred quite early (between 55Ma and 45Ma as I recall). If so the slab pull should have dropped very early. The previously subducting slab would still help drive mantle flow that could help entrain India though for a while.\nAuthor response: Thank you for your insightful comment, and we fully agree with your perspective. Indeed, the subducted continental crust/lithosphere is not a compelling candidate for the major driving force. Your point about the detached subducting slab driving mantle flow aligns with the main focus of our study. While our research does not directly disprove the continental slab pull theory, it demonstrates that mantle drag alone is sufficiently powerful to sustain the Cenozoic Tibetan motion without necessitating additional forces. Based on your comment, we have added the following statement in line 53-54: “ … , not to mention the likely early formation of eclogite densification18 and the early slab breakoff19.”\nTo provide further clarification, although our model takes into account the eclogization of the lower continental crust (new Supplementary Fig. 8), the dense eclogite does not substantially contribute to the formation of large-scale mantle flow. In our sensitivity analysis on the origin of lateral mantle flow, detailed in lines 109-126, all models have the same lithospheric structure, including the continental eclogite, while portions of the sinking slabs were removed to test their impacts on mantle flow. When all the slabs were absent (Supplementary Fig. 3), the rapid mantle flow disappeared, despite the presence of the lithospheric eclogite. This observation suggests that continental eclogites do not play a significant role in driving the mantle drag. We added in line 123-124: “The dense continental eclogite (see below) existing in all three cases have only minor impact on the large-scale mantle flow.” to further clarify this.\nLine 46. also at odd with geology as I just stated in my previous comment.\nAuthor response: Thank you very much for the comments. Following the previous comment, we updated this line (now line 53-54) with the geology evidences the reviewer mentioned.\nLine 50. What you mean here is not clear to me\nAuthor response: We want to emphasize that the strong force required to sustain the Tibetan Plateau is unlikely to be physically transmitted through the Indian-Capricorn plate boundary via shear stress. This is demonstrated in geodynamic models, such as the one presented in Zhou et al., 2020 (old ref. 20, new ref. 24). To better convey this point, we have modified line 57-58 to read, \"especially for a force large enough to overcome the substantial GPE-induced resistance from Tibet.\"\nLine 59. Now the pn you are addressing is much more clearly stated.\nAuthor response: Thank you.\nLine 65. Note that a recent article suggests that LLSVP would actually not be a slab graveyard but rather the remains of the impactor that for the moon (see article by Yuan Qian in the November issue of Nature). Does the interpretation of the origin of LLSVP actually matter? You don't necessarily need to discuss that possible new interpretation of LLSP. The proposed theory is very new and probably needs further investigation but i thought I would ask out of curiosity.\nAuthor response: Thank you for the interesting comment. In our model, the LLSVP is defined as a chemically distinct, uniform layer that is denser than the ambient mantle, and this definition is consistent throughout the model run (from 200 Ma to the present). In this sense, the conclusions of Yuan et al., 2023, do not directly contradict the structure used in our model.\nWe acknowledge that our starting point at 200 Ma might not be entirely realistic compared to the longevity of LLSVP. However, our model reproduces LLSVP shapes similar to those observed in seismic tomography. This similarity implies that the current shapes and locations of the LLSVPs may have primarily formed during the recent 200 Ma of mantle evolution and may not be strongly influenced by earlier events. Consequently, even considering the findings of Yuan et al., 2023, our conclusions remain valid within the context of our model setup and the timescale we investigated.\nLine 80. but not the slow down at the time of the collision. I assume this is because the resistive effect of Tibet's GPE is kept constant in your analysis, right?\nAuthor response: In our model, the GPE resistance is not artificially assigned but rather emerges as a self-sustained dynamic consequence of the Indian and Eurasian plates being pushed together by the driving forces. Prior to the collision, the Neo-Tethys Ocean separates the Indian and Eurasian continents, resulting in no collisional resistance being applied from the Eurasia side. The variation in resistance from the north of Indian plate over time is also depicted in the updated Fig. 4c, which demonstrates the model's ability to capture the evolving dynamics of the collision process.\nLine 89. I see a sublithospheric flow that is closer to Couette flow at 50Ma (with the Indian plate moving faster than the mantle) and a transition to a Poiseuille flow at 40Ma. India seems to have a very low surface velocity at this time and then it accelerates again to a larger velocity but still lower than upper mantle northward velocity so that a Poiseuille sub-lithospheric flow is maintained. Why is India motion so slow at 40Ma and why would it accelerate later on. This should not be happening in a NNR reference frame. What reference frame to you use to plot velocities in Fig 3. There is an inconsistency with the model velocity plotted in Fig 4 (yellow line).\nAuthor response: Thank you again for raising this important point, which we have addressed in a previous similar comment. To clarify again, Fig. 4c presents a combination of the assimilated plate reconstruction and the model output in the no-net-rotation (NNR) reference frame, ensuring a consistent comparison between the two datasets. On the other hand, Fig. 3d-f follows the reference frame used in the plate reconstruction by Muller et al., 2016 (old ref. 2, new ref. 3).\nLine 100. Good description\nAuthor response: Thank you.\nLine 105. I am puzzled by this statement. I follow the argument but velocities (so the plate velocities and upper mantle velocities) depend of the viscosity, so basal shear stress should depend quite a lot on the assumed viscosity structure.\nAuthor response: In the Stokes equation (Eqn. 2 in methods) we solve, the forces, including dynamic pressure gradients and stress, are conservative. Due to the given amount of gravitational force (determined by the density structure) present in the model, there is a fixed amount of pressure and local stress that can permit velocity changes. Although variations in viscosity may result in different velocity changes, the underlying force terms in the solution remain largely conservative. Put in a simplified mathematical form: in the mantle, the density anomaly of sinking slab creates the faster mantle flow following: velocity = density / viscosity, and close to the lithosphere-asthenosphere boundary we have mantle drag force following force = viscosity * velocity (more strictly its gradient), so force is only dependent on density, with viscosity largely canceled out. This is an important reason why we can estimate the magnitude of forces without worry too much about the velocity, as forces are conservative and derivable from density, while velocity is sensitive to the viscosity structure, which is less well constrained than density.\nLine 122. This is interesting, but I don't see why the slow down of India happens in your model. Also it looks like a very large slow down in Figure 3d and modest slow down in Figure 4...\nAuthor response: We apologize again for the confusion here. The apparent discrepancy in the magnitude of India's slowdown between Figure 3d and Figure 4 can be attributed to the different reference frames used, as mentioned in the previous section. In our model, the slowdown of the Indian plate from 50 Ma to 40 Ma is consistent with a similar decrease in the mantle drag during this period, providing a logical explanation for the observed deceleration. However, we acknowledge that this modeled slowdown does not fully match the reconstructed plate motion, as the rapid plate velocity at 50 Ma is not reproduced in our simulation. This suggests that additional factors or initial conditions not captured in our current model setup may have contributed to the high plate velocity at 50 Ma.\nLine 125. It seems to me that you are putting the cart before the horse. You can't claim that it is primary before discussing Figure 4.\nAuthor response: Thank you for highlighting this issue. To address your concern and improve the logical flow of the text, we have revised lines 133-136 as follows:” This fast asthenospheric flow persists from >50 Ma to the present day, suggesting a consistent driving force for the Indian plate throughout the post-collision Cenozoic era. In the following section, we quantify this driving force and demonstrate its dominant contribution to the Indian-Asian collision.”\nLine 138. I know that this terminology is common in the geodynamic literature but it bothers be a bit. In solid mechanics friction refer to shear stress proportional to normal stress along an interface.... What you call 'frictional stress' is the horizontal shear stress associated with sublithospheric shear flow. Not a big deal.\nAuthor response: Indeed, we use the term \"friction drag\" following conventions in the fluid dynamics community, which differs from the concept of \"frictional stress\" in solid mechanics. To avoid confusion, we have replaced the term \"frictional stress\" that you mentioned in a subsequent comment with \"shear stress,\" which more accurately describes the quantity accounted for. In solid mechanics, stress arises when shear strain occurs, and friction occurs when normal stress is applied, creating shear strain across contact boundaries. In fluid mechanics, friction occurs in the presence of a strain rate, which happens when fluid flows through contact boundaries. Although they differ in their constitutive context, the concepts of friction in solid and fluid mechanics bear similarities.\nLine 141. within the asthenosphere! (or below the lithosphere if you prefer). By (mechanical) definition of the lithosphere the vertical gradient of horizontal velocity is negligible within the lithosphere.\nAuthor response: Our analysis reveals that the maximum stress actually occurs at the base of the lithosphere, where the viscosity gradually transitions from strong to weak, yet remains stronger than the average asthenosphere. While this may not be the location of the maximum velocity gradient, the much higher viscosity in this region results in the maximum stress occurring within the stronger portion of the lithosphere. This finding is also mentioned in the updated line 173-175, where we state, \" … our results demonstrate that the maximum shear stress occurs within the relatively strong lithospheric root (Supplementary Fig. 4a)\"\nLine 164. This is not frictional! call it 'basal shear stress' may be.\nAuthor response: Thank you for bringing this to our attention. You are correct; it is not appropriate to describe the stress as frictional stress in this context. To address this, we have revised the text to use the term \"maximum shear stress\" instead in line 174.\nLine 165. To most people 'friction' implies a shear stress proportional to normal stress. The basal shear stress in your model is due to viscous shear flow.\nAuthor response: Thank you for emphasizing the importance of using appropriate terminology. As mentioned in our previous two responses, we agree that the term \"frictional stress\" is not suitable in this context, and we have revised the text accordingly. However, we would like to clarify that \"friction drag\" and \"pressure drag\" are well-defined concepts in the fluid dynamics context. Considering the time scale and constitutive fluid rheology of our simulation, we believe it is appropriate to describe the terminology within a fluid mechanics framework. Therefore, we think that retaining the term \"friction drag\" is justified in this case.\nLine 228. This doesn't seem so small (we got the same estimate in Copley et al) and of the same order of magnitude of the sa\nAuthor response: A ridge push force of 3×10^12 N/m, while not insignificant, is considerably smaller than the mantle drag force estimated from our model. Moreover, in the study by Ghosh et al., 2006, which provides the lowest estimate of Tibetan GPE (3×10^12 N/m), the ridge push force calculated using the same method is only 1×10^12 N/m, even smaller compared to the GPE. In both scenarios, ridge push is unlikely to be a major contributor to the India-Asia collision.\nLine 637. Why? I'd like to know the effect of slab break-off on the force balance.\nAuthor response: As discussed earlier, we have extended the resistance data presented in Fig. 4 to address your concern. In the figure caption, we have added the following statement to clarify the potential limitations of the northern resistance estimate before 40 Ma: “The purple curve shows the resistant force from continental collision in the north of the Indian plate. Note that it may not reflect full northern resistance before 40Ma given the slab is still attached to the lithosphere acting as a stress guide.” (line 768-770). This additional information acknowledges the potential influence of slab break-off on the force balance and the limitations of our current model in fully capturing its effect.\nFig. 4b You should report the forces and include the force transmitted and the boundaries with the Arabia and Capricorn plates.\nAuthor response: Thank you for your insightful comment. We have addressed this point in detail in the major comment section. Our current method has limitations in fully accounting for the force balance in the Arabian and Capricorn plates, as the force transmitted may be predominantly in the form of shear stress. That said, to investigate the potential influence of forces from the Capricorn plate, we have conducted additional tests, which indicate that these forces play a minor role in driving the collision. [Redacted] This finding suggests that the effects from the Arabian and Capricorn plates are not likely to be the primary drivers of the collision.\nReference\nZhong, S., Gurnis, M. & Moresi, L. Free-surface formulation of mantle convection-I. Basic theory and application to plumes. Geophys. J. Int. 127, 708–718 (1996).\nLi, Y., Liu, L., Peng, D., Dong, H. & Li, S. Evaluating tomotectonic plate reconstructions using geodynamic models with data assimilation, the case for North America. Earth-Sci. Rev. 244, 104518 (2023).\nCrameri, F., Lithgow‐Bertelloni, C. R. & Tackley, P. J. The dynamical control of subduction parameters on surface topography. Geochem. Geophys. Geosystems 18, 1661–1687 (2017).\nSimpson, R. W. Quantifying Anderson's fault types. Journal of Geophysical Research, 102(B8), 17909–17919 (1997).\nHeidbach, O., et al. The World Stress Map database release 2016: Crustal stress pattern across scales. Tectonophysics, 744, 484-498 (2018).\nCopley, A., Avouac, J. & Royer, J. India‐Asia collision and the Cenozoic slowdown of the Indian plate: Implications for the forces driving plate motions. J. Geophys. Res. Solid Earth 115, 2009JB006634 (2010).\nMüller, R. D. et al. Ocean Basin Evolution and Global-Scale Plate Reorganization Events Since Pangea Breakup. Annu. Rev. Earth Planet. Sci. 44, 107–138 (2016).\nBird, P. Testing hypotheses on plate‐driving mechanisms with global lithosphere models including topography, thermal structure, and faults. J. Geophys. Res. Solid Earth 103, 10115–10129 (1998).\nMolnar, P. & Lyon-Caen, H. Some simple physical aspects of the support, structure, and evolution of mountain belts. in Geological Society of America Special Papers vol. 218 179–208 (Geological Society of America, 1988).\nMolnar, P., England, P. & Martinod, J. Mantle dynamics, uplift of the Tibetan Plateau, and the Indian Monsoon. Rev. Geophys. 31, 357–396 (1993).\nForsyth, D. & Uyeda, S. On the Relative Importance of the Driving Forces of Plate Motion. Geophys. J. Int. 43, 163–200 (1975).\nZhou, X., Li, Z.-H., Gerya, T. V. & Stern, R. J. Lateral propagation–induced subduction initiation at passive continental margins controlled by preexisting lithospheric weakness. Sci. Adv. 6, eaaz1048 (2020).\nYuan, Q., et al., Moon-forming impactor as a source of Earth’s basal mantle anomalies. Nature, 623(7985), 95-99 (2023).\nGhosh, A., Holt, W. E., Flesch, L. M. & Haines, A. J. Gravitational potential energy of the Tibetan Plateau and the forces driving the Indian plate. Geology 34, 321 (2006).\nReviewer #2\nThis is an interesting study in which the authors use a classical geodynamic mantle convection code, CitcomS, to make the case that large scale mantle flow is what keeps the Tibetan Plateau up.\nThe idea that a large scale mantle wind drives plates is not new; Faccenna and Becker (2010 Nature) have popularized this for the Mediterranean and applied it later to the India-Asia collision zone in Becker and Faccenna 2011 (using the same code). The results in the current study largely confirm this and also demonstrate that slabs in the lower mantle play a key role in driving the plates. The study is an improvement over previous work because 1) a novel data assimilation method is employed that takes plate reconstructions into account; 2) a more serious effort is done to quantify the forces that occur in the simulations.\nThis is clearly of interest to a broad audience, and the codes and data assimilation method have been used in many previous studies.\nAuthor response: Thank you very much for the encouraging comments.\nYet, I do have two major concerns, that can be addressed with additional simulations, as well as a few more minor remarks:\n1) A key aspect of this paper is that the gravitational potential energy of the Tibetan plateau, which is caused by the difference in elevation of the plateau and the surrounding lower areas, is very significant and on the order of 3-9x1012 N/m as discussed in the main text. The goal of the paper is to quantify mechanisms that are able to provide sufficient forces and the authors do a good job in that. Yet, a key limitation of their models is that it does not include a free surface and as a result, the gravitational potential energy (the resistive force) is not taken into account. The authors are aware of this limitation and point out that a free surface is not being used (lines 160-163). Yet, they only argue on the basis of the forces that are estimated from natural observations, while using the driving forces from their model without topography. In models with an actual free surface, the magnitude of these driving forces may actually be quite different from those of free slip models; I am not aware of a study that has quantified this error and therefore I cannot judge whether the key conclusions of this study will hold. If you implement a setup with a changing crustal thickness and a free slip top boundary, the crust wants to move upwards due to a density instability that causes additional forces. As such\nWhereas I realize that it is unfair to ask the authors to change CitcomS to take a real free surface into account, one could estimate the error by performing 2D or 3D idealized simulations with other open-source geodynamic codes (such as ASPECT, pTatin3D, or LaMEM) that do have a self-consistent free surface implemented. It would be rather straightforward to do a similar force balance analysis for a setup with a craton that has roots of varying thickness and a large scale mantle flow, and uses either free slip or free surface upper boundary conditions. Doing\nAuthor response: Thank you for bringing up this crucial point regarding the potential limitations of using a free-slip surface boundary condition in CitcomS. We appreciate your insight, as the other reviewer has also expressed similar concerns. To address this issue comprehensively, we would like to provide not only a detailed explanation of our reasoning but also the results of a preliminary ASPECT model test that we have conducted.\nWe agree that a free-slip model has some difference from a free surface model in capturing the gravitational potential energy (GPE). However, we believe that our key findings regarding the plate-wise force balance, particularly the magnitude of mantle drag and the resistance related to GPE, remain robust for the following reasons:\n1. Previous studies have shown that the difference in mantle flow between free-slip and free-surface models is relatively minor (Zhong et al., 1996). This suggests that the resulting plate-driving forces, especially due to mantle drag, should be similar between the two types of models. Mechanically, a free-slip boundary, although with a flat surface, stores a finite vertical stress component that is equivalent to the stress as the same depth of a free-slip model; this means their respective GPE distribution is also similar. Indeed, codes like CitcomS compute surface topography by converting this sub-surface stress into a topography that very much resemble the real case, as seen from the figure below (Review Fig. 4) that we just published.\nReview Figure 4. Predicted topographic evolution of North America based on a data-assimilation model from CitcomS (Li et al., ESR, 2023). The modeled topography resembles the real Earth topography both in the amplitude (>5 km) and geological observations (red arrows at 120 Ma mark paleo current direction, and red lines at 80 Ma represent paleo shorelines).\n2. Some earlier studies show that free-slip models can generate more short-wavelength oscillations in topography compared to free-surface models (e.g., Fig. 4 in Crameri et al., 2017), but the longer-wavelength topographic features are generally consistent. Zhong et al. (1996) show that the smoothing in topography in free-surface case is due to the temporal retard from surface relaxation. In the Crameri et al. (2017) study, the slab geometry is somewhat artificially defined, i.e., not naturally produced by the model. We suspect, if the slab was generated more naturally (thus better trained by the boundary condition), the calculated (from sub-surface stress) surface topography would be more consistent with the free-surface case, as shown below in our new test. In fact, even assuming such short-wavelength topographic differences do exist, they would have a small effect on the overall plate force distribution. Since our analysis focuses on the plate-scale force balance rather than local topographic variations, these short-wavelength differences are unlikely to significantly impact our overall conclusions. On top of this, such short-wavelength oscillation is also related to the rough surface viscosity variation in similar spatial scale (e.g., Fig. 4a,b,d,e in Crameri et al., 2017, likely due to plastic deformation that amplifies the intrinsic rheological inconsistency with the arbitrarily defined slab geometry), which is not the case in our model (Fig. 3c).\n3. To further test the robustness of our results, we conducted an idealized numerical experiment using the geodynamic modeling software ASPECT. We ran both free-slip and free-surface models of mantle-drag-driven continental collision over 0.4 Myr. As expected, the free-surface model produced a smoother topographic profile while the free-slip model exhibited more short-wavelength features (Review Fig. 5b). Crucially, however, the vertically integrated lithospheric pressure, which represents the accumulation of driving and resisting forces, was very similar between the two models (Review Fig. 5a). Additionally, the mantle flow patterns and plate velocities were nearly identical (Review Fig. 5c, d). All these confirm our previous inference and support the force analysis in a free-slip model.\nThese results give us confidence that, although our CitcomS models do not explicitly include a free surface, they still capture the first-order effects of GPE on the lithospheric stress state and plate-driving forces. The strong similarity in vertically integrated pressure between the free-slip and free-surface ASPECT models supports the validity of the resistive force we infer along the Himalaya plate boundary. We have added in line 319-336 “Previous research suggests that mantle models based on a free-slip… as it captures the essential features of the stress distribution” to elaborate on this in the main text, and included Review Fig. 5 as our new Supplementary Fig. 10.\nWe acknowledge that future work is needed to fully establish the force variation between a free-slip and free-surface model being minor. In this sense, the additional ASPECT models are beyond the scope of the current manuscript and will be developed into a separate publication.\nReview Figure. 5 parison of results between the free slip and free surface models. (a) Depth integrated lithospheric dynamic pressure and (b) topography for both the free-slip and free-surface models. Lithospheric dynamic pressure is integrated vertically from the surface to 1000°C isotherm within the lithosphere, which is roughly equivalent to the 1e21 viscosity boundary. It directly reflects the lithospheric force accumulation along the horizontal direction. (c) Free-slip and (d) free-surface models showing viscosity and velocity.\n2) The resolution of the models is rather coarse with only 113 nodes in the vertical direction, and a vertical resolution of around 12 km in the crust. This is hardly 4 elements over the whole crust which is on the lower bound of what is required to even resolve this (given the linear elements that CitcomS uses). It is clearly crucial that the authors perform resolution sensitivity tests to demonstrate that their results are robust.\nAuthor response: This is a good point. While it is true that we have a limited number of elements representing the crust due to the coarse vertical resolution, it is important to note that the stress accumulation in our model occurs primarily within the lithosphere as a whole rather than the crust alone. Although we use 5 separate compositions (Supplementary Fig. 1d, Supplementary Table 3) to define the continental crust and continental mantle lithosphere, the viscosity cut-off results in the entire continental lithosphere interior having mostly uniform viscosity (see Fig. 3c, for example), making it mechanically coherent. Since our entire continental lithosphere consists of 17 nodes vertically, it is sufficient to capture the essential large-scale force balance.\nWe would also like to further clarify about the horizontal resolution used in our model, which we believe is also able to capture essential large-scale force balance. This is because of the nature of assimilated rigid-plate velocity that prohibits severe lithospheric shortening within the interior of Tibet. Therefore, the deformation is more localized near the plate boundary. That said, the force balance across the convergent boundary is realized properly as the driver is from the southern side, and this is why we evaluate the plateau resistance within the northern edge of the Indian plate.\nOn top of this, recent works from our group (e.g., Cao and Liu, 2021, as well as the new ref. 29-33) utilize models with similar spatial resolution for Western North America, which has an thinner lithosphere compared to that in India and Tibet, but still capture the essence of 3D stress and reproduce topography and crustal deformation. This further confirms the validity of the resolution we used. We added further elaboration in line 411-417 “The entire continental lithosphere including the crust and the mantle lithosphere within a plate is mechanically coherent because of viscosity cut-off (Fig. 3c). This allows the stress accumulation in our model occurs primarily within the lithosphere as a whole … recent studies29-34 utilizing similar model resolution successfully capture 3D lithospheric stress state, deformation and topography.” We hope this addresses the reviewer’s concern.\nMinor remarks:\n- For how long have the forces been sufficient large to keep the Tibetan plateau high?\nAuthor response: Based on the results presented in Fig. 4c, the mantle drag force has been sufficiently high since 50 Ma, and the resistance force directed towards the northern boundary has been sufficiently high since 40 Ma, shortly after the collision. These findings suggest that the forces have been large enough to sustain the elevated Tibetan Plateau since the early stage of the India-Eurasia collision.\n- Please use SI units throughout (2800 kg/m3 instead of 2.8 g/cm3)\nAuthor response: Thank you for pointing out the inconsistency in units. We have corrected the density unit to \"2800 kg/m3\" in line 407 to adhere to SI units. However, we would like to note that units such as \"cm/yr\" and \"km\" are widely used and accepted in solid Earth studies. We believe that converting these units to their SI equivalents (e.g., \"m/s\" and \"m\") might not necessarily improve the overall readability of the manuscript given their scale. Nonetheless, we have carefully reviewed the manuscript to ensure that all other units are consistently presented and adhere to SI units where appropriate.\n- Please convert viscosity values on Fig. 3c into dimensional units (Pas)\nAuthor response: Thank you for your suggestion. We have updated Fig. 3c to present the viscosity values in dimensional units (Pa s) to ensure consistency and clarity throughout the manuscript.\n- Please add a table with all material parameters employed. The statement “more model information is available in recent studies using similar methods” is very unsatisfying as it doesn’t allow me to reproduce what was actually done.\nAuthor response: We apologize for the inconvenience caused by the lack of detailed information on the material parameters used in our study. To address this issue and to ensure reproducibility, we have now included three supplementary tables (Supplementary Tables 1, 2, and 3) that provide a comprehensive list of all material parameters employed in our models. These tables cover the relevant parameters for each component of the model, allowing for a clear understanding of the input values and enabling other researchers to reproduce our results. We appreciate your feedback, as it has helped us improve the clarity and completeness of our manuscript. Please refer to the newly added supplementary tables for the specific material parameters used in our study.\nReference\nZhong, S., Gurnis, M. & Moresi, L. Free-surface formulation of mantle convection-I. Basic theory and application to plumes. Geophys. J. Int. 127, 708–718 (1996).\nLi, Y., Liu, L., Peng, D., Dong, H. & Li, S. Evaluating tomotectonic plate reconstructions using geodynamic models with data assimilation, the case for North America. Earth-Sci. Rev. 244, 104518 (2023).\nCrameri, F., Lithgow‐Bertelloni, C. R. & Tackley, P. J. The dynamical control of subduction parameters on surface topography. Geochem. Geophys. Geosystems 18, 1661–1687 (2017).\nCao, Z. & Liu, L. Origin of Three‐Dimensional Crustal Stress Over the Conterminous United States. J. Geophys. Res. Solid Earth 126, e2021JB022137 (2021)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I have read the revision and the rebuttal letter. The additional test comparing simulations with a free-surface and free-slip boundary show that the two boundary conditions are, not surprisingly, not equivalent but that the impact of the difference on the simulations is small. I am now more convinced that the results presented in the study are solid. I am also satisfied by how the authors have addressed my other concerns, which were more minor. The manuscript is well written and the places where the model is not completely satisfying are discussed honestly. It makes a good contribution.\nI am pleased to recommend that the revised mansucript be accepted for publication.", "reviewer_2": "I am overall quite satisfied with the revisions, particularly with the additional simulations done to compare free slip vs. free surface simulations. A few minor remarks remain:\n1. In arguing for free slip vs. free surface, you actually simulate a different geometry than what is shown on for example figure 3c of the main text, where no weak zones in the crust are resolved (if weak zones are not resolved slab pull is not efficiently transmitted so one may artificially increase the relative importance of basal drag).\n2. In addition to the point above, to really make your case you should also include CitcomS models at the same resolution as your main simulations in the comparison plot shown in supplementary figure 10."}, "author_response": "We thank both reviewers for the kind support and reviewer #2 for the constructive follow up comments. These comments make good points and we have responded to each of the comments one-to-one. The original comments are in black font and our responses in blue.\nWe hope the additional work addresses all the reviewer’s concerns. The following are detailed one-to-one responses.\nReviewer #1\nI have read the revision and the rebuttal letter. The additional test comparing simulations with a free-surface and free-slip boundary show that the two boundary conditions are, not surprisingly, not equivalent but that the impact of the difference on the simulations is small. I am now more convinced that the results presented in the study are solid. I am also satisfied by how the authors have addressed my other concerns, which were more minor. The manuscript is well written and the places where the model is not completely satisfying are discussed honestly. It makes a good contribution.\nI am pleased to recommend that the revised manuscript be accepted for publication.\nAuthor response: Thank you very much for your encouraging comments and support.\nReviewer #2\nI am overall quite satisfied with the revisions, particularly with the additional simulations done to compare free slip vs. free surface simulations.\nAuthor response: Thank you very much for your encouraging comments and support.\nA few minor remarks remain:\n1. In arguing for free slip vs. free surface, you actually simulate a different geometry than what is shown on for example figure 3c of the main text, where no weak zones in the crust are resolved (if weak zones are not resolved slab pull is not efficiently transmitted so one may artificially increase the relative importance of basal drag).\nAuthor response: Thank you for bringing up this question. We agree that if weak zones are not resolved, the amount of mantle drag may artificially increase. However, we believe that in Fig. 3c, the lithospheric weak zone between the India and Eurasia plates is sufficiently weakened, including both the crustal part and the lithospheric mantle part. Indeed, the viscosity of the crustal part of the weak zone (1020 Pa s) is equivalent to the upper mantle viscosity in the model, ensuring effective decoupling. To further demonstrate that the weak zone is sufficiently weakened, we have included a von Mises stress cross-section (Review Fig. 1) which shows that there is nearly no deviatoric stress within the entire lithospheric weak zone. Had the crustal part not been weakened, there will be stress accumulation just like that in the plate interior, and there will be less realistic plate motion differentiation. We have added “, while the weak zone fully decouples the plate boundary between two plates” in line 415 and “Note that the collisional plate boundary between India and Eurasia at ~30˚ N is fully weakened in the entire lithospheric depth.” in Fig.3 caption, line 757-758 to clarify this.\nReview Figure 1. Deviatoric stress distribution of India-Asia collision system. The cross section corresponds to the location of Fig. 3b and 3c at 0 Ma. The blue arrows represent the modeled deviatoric principal stress, while the dark contour indicates the 0.6 non-dimensional temperature volume. At 30° N, the absence of von Mises Stress throughout the entire lithospheric depth suggests a completely weakened plate boundary between India and Eurasia.\nFurthermore, since there is no actively sinking slab mechanically attached to the lithosphere as shown in Fig. 3c, there is no slab pull competing with mantle drag for the Indian plate motion. This is also illustrated with our discussion of Neo-Tethyan slab detachment on the northern Indian plate boundary. Therefore, in this specific case, the relative importance of mantle drag to slab pull is not artificially increased.\nWe acknowledge that the weak zone in the ASPECT model has an inclined geometry, while the CitcomS model has a vertical boundary, which is a simplification under the coarser resolution. However, the CitcomS model's weak zone is 100 km wide, as compared to 12 km in the ASPECT model. This can ensure the two plates are fully decoupled in the CitcomS model. Thus, we believe that our CitcomS model has a sufficient weak zone to accounf for the force balance.\n2. In addition to the point above, to really make your case you should also include CitcomS models at the same resolution as your main simulations in the comparison plot shown in supplementary figure 10.\nAuthor response: The comparison you mentioned is indeed useful. However, the ASPECT model (Supplementary Fig. 10) can reach its current resolution due to two key factors: 1) it is a two-dimensional model, and 2) it employs adaptive mesh refinement. These features are not available in our CitcomS model, which is a three-dimensional global model without adaptive mesh refinement. Incorporating these features in the CitcomS model would require significant computational resources. To put this into perspective, the current ASPECT model takes around half a day to run, and each snapshot of the free-slip version for our global CitcomS model also takes a similar amount of time, and with 5376 CPUs! Given our current computational capabilities, it is not feasible for the CitcomS model to reach the same resolution as the ASPECT model. In fact, our global already has the highest published lateral resolution (25 km) among its peers.\nWhile it is possible to run a lower resolution ASPECT model to closely approximate the CitcomS model situation, we believe that it is unnecessary for the scope of this paper. Firstly, the large-scale topography and forces (Supplementary Fig. 10) occur at scales much larger than the horizontal resolution in the CitcomS global model (25km). This means that all the critical large-scale features demonstrated in the ASPECT models should be well reproduced by the CitcomS model at its current resolution. Secondly, the purpose of running the ASPECT model comparison is not to completely reproduce the CitcomS model result but to demonstrate that the force distribution within the lithosphere is similar between free-slip and free-surface models. For this purpose, we think the current model is sufficient. A more in-depth, thorough comparison between the two types of boundary conditions would be better suited for a future paper, which we plan to write. To acknowledge the reviewer’s suggestion, we have added “The overall similarities between topography and forces are of a wavelength much larger than the resolution of the global model (25 km).” in line 331-332 and “Based on these findings, we conclude that our global free-slip model is…” in line 336 to further clarify this idea."}]

Decision letter