Impact of Decelerating India‐Asia Convergence on the Crustal Flow Kinematics in Tibet: An Insight From Scaled Laboratory Modeling

The factors controlling the spatiotemporally varying deformation patterns in Tibet, a prolonged period (∼50 to 19 ± 3 Ma) of NNE‐SSW shortening, accompanied by eastward flow and orogen‐parallel extension in a later stage (19 ± 3 Ma to present‐day), are still poorly constrained. Using viscous models, we performed scaled laboratory experiments with steady and unsteady state collision kinematics to address this issue. Our model Tibet under steady‐state collision, irrespective of high (5.5 cm/yr) or low (3.5 cm/yr) indentation rates fails to produce the present‐day crustal velocity fields and the deformation patterns, reported from GPS observations. An unsteady‐state collision with decelerating convergence rates (5.5–3.5 cm/yr) is found to be a necessary condition for the initiation of eastward flow and ESE‐WNW extensional deformations. The model results also suggest that the mechanical resistance offered by the rigid Tarim block resulted in crustal uplift at faster rates in western Tibet, setting a west to east topographic gradient, existing till present‐day. This topographic gradient eventually polarized the gravity‐controlled flow in the east direction when the convergence velocity decelerated to ∼3.5 cm/yr at around 19 ± 3 Ma. Our model shows the present‐day eastward flow in central Tibet follows nearly a Poiseuille type velocity profile, bounded by the Himalaya in the south and the Tarim basin in northern Tibet. This flow kinematics allows us to explain the preferential locations of crustal‐scale dextral and sinistral faults in southern and northern Tibet, respectively. Finally, the present‐day model crustal‐flow velocity, strain‐rates, and topographic variations are validated with GPS and geological field data.


10.1029/2021GC009967
3 of 26 of the plateau and steep topography in its eastern and southern margin (Bird, 1991;Clark & Royden, 2000;Zhao & Morgan, 1987), and (f) an outward growth model from a pre-existing proto-Tibetan plateau (Murphy et al., 1997;. Yin (2000) proposed a completely different geodynamic model to explain E-W extension in the late Cenozoic time. According to his model, slab rollback of the western pacific plate gave rise to a regional lithospheric extension in entire southeast Asia, leading to Tibetan and Baikal rifting, and Shanxi graben formation. Despite a strong divergence in the prevailing views, it is generally accepted, based on the spectacular records of Tibet's long-term evolution history, that the evolution of plateau topography had an intricate relationship with the continental scale deformation episodes in Tibetan crust (DeCelles et al., 2002;England & Houseman, 1989; Y. Li, Wang, et al., 2015). Recent tectonic models have shown that the lateral crustal heterogeneities, such as rigid Tarim and Sichuan basins (Figure 1a) in the Asian plate have strongly influenced the topographic uplift and the long-term crustal deformations in the Tibetan plateau (L. Chen et al., 2017;Pusok & Kaus, 2015;Y. Yang & Liu, 2013).
Paleomagnetic data and plate reconstructions suggest that the India-Asia collision slowed the convergence rates, from 15 cm/yr at 50 Ma to 3.5 cm/yr at the present-day (Figure 1d) (Copley et al., 2010;Molnar & Stock, 2009;Pusok & Stegman, 2020). Such a decelerating convergence motion had the potential to reduce the tectonic forcing upon the deforming Asian plate, resulting in the gravitational collapse of the elevated plateau in the course of collision Rey et al., 2001). Copley and McKenzie's theoretical  (Gan et al., 2007). (b) Topographic profiles along CD (east-west), EF (northwest-southeast), and GH (southwest-northeast) sections (profiles are marked with dotted lines in panel a), showing a general west-to-east topographic gradient in the presentday Tibetan plateau. (c) Orogen-transverse variation of the east-west (N110°E) GPS velocity components, measured along N-S profiles (marked in solid colors in panel a). (d) Variation of the India-Asia convergence velocity with time. Note that different stages (S1, S2, and S3) with corresponding indentation velocity (cm/ yr) and crustal shortening (ϵ) are for the Type 3 experiment (discussed later). model (2007) predicts this gravitational collapse as the driving factor in setting the present-day extensional tectonics in Tibet. On the other hand, several workers suggested eastward gravitational spreading of the whole Tibetan lithosphere, confined between two rigid blocks (Indian craton in the south and the Tarim and Qaidam basins in the north) controlled the rifting and strike slip faulting in the Tibetan plateau (Bajolet et al., 2015;Fournier et al., 2004;Yin & Taylor, 2011). From laboratory models, Bajolet et al. (2013) showed that the symmetric lateral extrusion between two rigid colliding plates, coupled with a late stage gravity collapse of the thickened crust modulated Tibetan tectonics, including the syntaxis formation on the eastern and western flanks of Tibet. This model is in agreement with Fournier et al.,'s (2004) study, suggesting a free or weakly confined eastward boundary of Tibet as a necessary condition for eastward gravitational spreading of the Tibetan crust. However, these earlier numerical and analogue experiments mostly ignored the temporally varying India-Asia collisional kinematics in their models, setting the experimental runs at a constant boundary condition, either with an average convergence velocity of 5 cm/yr (Beaumont et al., 2001;Y. Yang & Liu, 2013) or an average indentation rate of ∼3.5 cm/yr (Schellart et al., 2019) for the entire 50 Ma collision event. Furthermore, most of the earlier laboratory models did not account for lateral crustal heterogeneities, particularly rigid Tarim and Sichuan basins situated in the north-western and easternmost regions of Tibet, respectively. Numerical simulations of England and Houseman (1989) took into account the effects of decreasing collision force, expressed in terms of Argand number, but their study excludes any possible effects of the lateral crustal heterogeneities (Tarim and Sichuan basins) around the Tibetan plateau. Despite a remarkable progress of tectonic modeling in recent time, it is yet to address the following questions-did the temporally varying collision kinematics play a critical role in the tectonic evolution of Tibet, and if so, what were the complementary effects of the lateral crustal heterogeneities? These issues constitute the central theme of our present article.
From indentation model experiments, we hypothesize that the decelerating convergence motion in India-Asia collision and the presence of lateral crustal heterogeneities are the most crucial factors in switching the tectonic transition in the Tibetan plateau, an initial period (∼50 to 19 ± 3 Ma) of NNE-SSW crustal shortening followed by initiation of eastward crustal flow and ESE-WNW extensional tectonics (19 ± 3 to 0 Ma), as observed in its present-day setting. To test this hypothesis, we considered three kinematic states: (a) constant fast-rate collision, (b) constant slow-rate collision, and (c) decelerating collision rate in running the laboratory model experiments. Our model experiment with the decelerating collision kinematics reproduces the current eastward crustal flow pattern in central Tibet revealed from the GPS studies. The N-S profile of the present-day eastward flow components is shown to follow Poiseuille behavior, where the Himalaya in the south and the Tarim basin in the north act as the guiding walls. Finally, we compare our model results with the present-day surface deformation pattern, estimated from GPS data and geological records in the Tibetan plateau, and provide a new tectonic evolution model for the plateau and its surrounding regions.

Modeling Approach
We adopted a classical indentation tectonic model to perform analogue experiments of the India-Asia collision in laboratory (Fournier et al., 2004;Schellart et al., 2019;Tapponnier et al., 1982). The experiments were designed to study the Tibetan plateau deformation in the last ∼50 Ma as a function of the three collisional kinematic states discussed in the preceding section. We use the experimental results to constrain the decelerating India-Asia collision rate as a critical collision kinematics required for setting the eastward crustal flow in Tibet. Our laboratory models allow us to explain the crustal velocity pattern, topography, and strain fields in space and time.

Model Setup
The continental lithosphere is often modeled as a fluid continuum of Newtonian or power-law rheology to describe its large-scale deformation behavior on million years' time scales (England & McKenzie, 1982). As the horizontal length scale of the continuum system (Tibetan plateau) far exceeds the lithospheric thickness, several workers have developed their tectonic models within a framework of thin viscous sheet (TVS) approximation (Bischoff & Flesch, 2018;Yang & Liu, 2013). The TVS approximation (England & McKenzie, 1982) simplifies depth-wise rheological variations in the lithosphere by considering a vertically averaged and approximated lithospheric rheology, and treats the system as a single continuum medium. We chose a two-layer continuum model to perform scaled laboratory experiments on Tibet.
In our model, the continental deformation in the Asian domain has been enforced by imposing external velocities in the Indian indenter, attached to a step motor ( Figure 2a). It is noteworthy that the dynamics of convergent margins involves a balance of the far-field tectonic forces, orogenic loading and other forces (Clark, 2012;Iaffaldano et al., 2006;Regard et al., 2005), which eventually determines the kinematics in orogens. However, to investigate exclusively the effects of reducing convergence rate on the mode of Tibetan crustal deformations, we chose a kinematic model boundary condition by externally controlling the velocities at the model boundaries. This kinematic approach to modeling is similar to those used in many earlier India-Asia collision models (Bajolet et al., 2015;Bischoff & Flesch, 2018;Cook & Royden, 2008;Fournier et al., 2004;Maiti & Mandal, 2021;Schellart et al., 2019;Y. Yang & Liu, 2013). We simulate approximately north-east ward Indian indentation into Asia. The northern boundary (leading) of the Indian indenter represents the contact between the Indian plate and the Himalaya-Tibetan plateau. A recent analysis of the seismic tomography and earthquake locations data suggests that the underthrusting of the Indian plate beneath the Tibetan plateau reaches up to ∼300-500 km north of the Himalayan front (Craig et al., 2020). To include this subduction induced simple shear beneath the Asian plate, we have designed the Indian indenter, allowing a part of it to underthrust for about 350 km north of the mountain front ( Figure 2b).
Unlike previous laboratory experiments of India-Asia collision (Bajolet et al., 2015;Fournier et al., 2004;Schellart et al., 2019;Tapponnier et al., 1982), we have incorporated several important additional features in our present model. One of them replicates the presence of lateral crustal heterogeneities that is, the inclusion of rigid Tarim and Sichuan basins, situated in the northern and the eastern margins of the Tibetan plateau, respectively. They act as a region of strong lithosphere (Clark et al., 2005;Dewey et al., 1988; Penney & Copley, 2021), as predicted from observations of fast seismic velocities down to a depth of 200 km (Li et al., 2006), measurements of low heat flow rate (Wang, 2001), and present-day GPS velocity data (Gan et al., 2007;Liang et al., 2013;Zhang et al., 2004). On the other hand, we have considered an arcuate plate boundary of the Indian plate. It is noteworthy that all the previous analogue India-Asia collision models used the Indian indenter ideally rectangular in shape, forming a straight plate contact. But, the present-day southern margin of the Tibetan plateau, that is, India-Asia contact, is strongly curved. The plate reconstructions and rotations measured from palaeomagnetic data (as summarized in Van Hinsbergen et al., 2011) suggest that this arcuate geometry existed since the onset of continent-continent collision in this region. Copley (2012) has shown how this curved plate contact influenced the gravitational spreading of the Tibetan plateau in the late stage of the orogeny. As our present study primarily aims to investigate the entire tectonic evolution of the plateau, including its late-stage gravity collapse, we choose such a curved plate contact to model the India-Asia collision system.
There are a number of advantages in analogue modeling used in this study. (a) Analogue models provide us an excellent scope for studying the four-dimensional (three-dimensional space + time) evolution of Tibet. (b) They have a free top surface, allowing the continuum to deform and produce mountain topography and plateau morphology, which can be easily compared with the natural prototype that is, the Himalaya-Tibet Mountain system. (c) The scaled analogue models reproduce large scale (i.e., 3,500 km × 5,000 km continental lithosphere) crustal flows, which can be studied in high resolution to compare the model results with natural data.

Material Properties and Scaling
We have scaled our analogue model to the natural prototype setting, keeping geometric similarity of the model domain, similarity in the temporal variation of collision kinematics, and downscaling of the crustal and mantle rheology (Hubbert, 1937;Ribe & Davaille, 2013;Schellart & Strak, 2016). To model the India-Asia collision system we have used a length scale factor (nature to model ratio): l r = 1.25 × 10 −7 . This scaling yields a model dimension of 40 cm to represent an approximately 3,200 km east-west length of the Tibetan plateau. Similarly, a 40 km thick crustal lithosphere in the initial model is represented as 0.5 cm thick viscous layer. The model crustal lithosphere rests over an initially 2 cm thick layer, which represents a sub-crustal mantle region. The whole setup is contained within an experimental tank. At the contact between the model base and the basal plate of the experimental tank, we have used liquid soap (viscosity ∼0.5 Pa s) to minimize friction to the model. To ensure the lubrication at the desired level, we performed independent experiments and evaluated base normal velocity profiles, calculated from particle image velocimetry (PIV) on the model side faces (details provided in Figure S2 in Supporting Information S1).
Model materials with different rheological properties were chosen to appropriately scale the crustal lithosphere and underlying sub-crustal mantle region in the experiments. A glucose syrup (syrup-1) of constant density (ρ cl = 1,406 ± 3 kg/m 3 ) and viscosity (µ cl = 305 ± 5 Pa s) is used to model the whole (40 km thick) crustal-lithosphere as a single continuum, representing its depth-averaged physical properties in the corresponding natural prototype (viscosity: 3 × 10 21 Pa s; density: 2,780 kg/m 3 ) (Jiménez-Munt et al., 2008). The model to nature scaling factors for viscosity (µ r ) and density (ρ r ) are found to be 10 −19 and ∼0.5, respectively. We have used another variety of glucose syrup (syrup-2) with different density and viscosity (µ sm = 28 ± 5 Pa s; ρ sm = 1,510 ± 5 kg/m 3 ) to simulate the sub-crustal mantle region (viscosity: ∼10 20 Pa s; density: 3,200 kg/m 3 ) (Jiménez-Munt et al., 2008). To obtain the right density of glucose syrup-2, we have mixed water soluble sodium carbonate (Na 2 CO 3 . 10 H 2 O, ρ = 2.53 gm/cm 3 ) with a low-viscosity glucose syrup (µ = 10 ± 3 Pa s and ρ = 1,392 ± 5 kg/m 3 ). The lower µ sm and higher ρ sm of syrup-2 than that of syrup-1 allows the base of the crustal lithosphere to deform vertically and maintain an isostatic equilibrium of the model topography during the experimental run. We have chosen commercial plasticine mixed with oil (viscosity = 3 × 10 5 Pa s; density = 1.5 gm/cc) (Zulauf & Zulauf, 2004), which is relatively much stronger than the glucose syrup-1, to represent rigid Tarim and Sichuan basins in the model. The viscosity ratio in the model is then in the order of ∼10 3 , which agrees well with the numerical estimates of the viscosity ratio (∼10 2 to 10 3 ) between rigid Tarim basin and Tibetan upper crust (Clark & Royden, 2000;Flesch et al., 2001;Penney & Copley, 2021).
Considering an extremely slow motion in the geodynamic setting, the dynamic scaling must satisfy the following relation (Ghosh et al., 2020;Marques & Mandal, 2016), Equation 1 yields a time ratio: t r = 1.6 × 10 −12 . That means, 1 Ma time in nature equals to ∼50 s in the scaled model. Now, using the relation: the velocity ratio in our model is: U r = 7.813 × 10 4 . According to this scaling factor, Indian plate's indentation velocity of 5.5, 4.5, and 3.5 cm/yr in nature equates to 8.17, 6.69, and 5.20 mm/min in the model, respectively. The total model runtime in our experiments was ∼42 min, which is equivalent to ∼50 Ma in nature. All the scaling parameters are listed in Table 1.

Recording of Experiments
For a quantitative analysis of the model crustal flow velocity and topography development in each experiment, we used a particle image velocimetry (PIV) system, which consisted of four synchronized cameras systematically aligned around the experimental box: three at the top and one at the side of the model set-up ( Figure S3 in Supporting Information S1). The top cameras were placed ∼80 cm above the set-up table, with 10°-15° angles between them to capture sharp images (Galland et al., 2016). All four cameras were standard Digital single-lens reflex cameras (Nikon D3400 24.2 MP) with a shutter speed of 1/125 and a focal length of 24 mm. The operation of these four cameras was synchronized with the help of open-source software, digiCamControl in a single computer system.
We used the structure-from-motion technique to retrieve quantitative surface deformation in the experiments (Galland et al., 2016). A stand-alone spatial photogrammetry software Micmac (Galland et al., 2016;Girod, 2012) was used to process the photographs and create a series of digital elevation models (DEMs) of the evolving topographic surface. The high-resolution DEMs were constructed after defining ground control points and the scale between the model and the images (Galland et al., 2016). The DEMs from the successive snaps of a model run gave a detailed time-series configuration of the evolving surface topography. We adopted this approach for a quantitative analysis of the topographic evolution in the Tibetan Plateau.
The method of high-resolution digital image correlation was implemented to track and analyze the evolving displacement, stress and strain fields in our scaled analogue models. Some raw photographic images of our model experiments are provided as examples in supporting information ( Figure S4 in Supporting Information S1). We used PIVLab (Thielicke & Stamhuis, 2010) Note. The model scaling factors are provided in the extreme right column.  (Cook & Royden, 2008;Jiménez-Munt et al., 2008;Liu & Yang, 2003;Schellart et al., 2019) images are adjusted for topographic relief, lens distortion and camera tilt and they provide an accurate representation of the photographed surface. The stress and strain fields were deduced with high accuracy. We performed several such quantitative analyses of the model observations to study the effect of the India-Asia convergence kinematics on the Tibetan Plateau.

Reference Experiments
Plate reconstructions studies suggest that at the time of India-Asia collision, ∼50 Ma ago, the Indian subcontinent was located at ∼2,500 km south of its current position (Johnson, 2002;Molnar et al., 1993;Rowley, 1998;Royden et al., 2008;Van Hinsbergen et al., 2011). The amount of shortening involved in the evolution of the Himalaya-Tibet system is, however, a complicated issue. Again, there is a serious uncertainty with the partitioning of total shortening into the Himalayan range and Asia as the available estimates of Asian shortening vary on a wide spectrum (Ingalls et al., 2016;Johnson, 2002;Van Hinsbergen et al., 2011;Yakovlev, 2015). Recently, Schellart et al. (2019) demonstrated from experimental studies that 1,740 ± 300 km of India-Eurasia convergence since the beginning of collision has been accommodated by indentation driven Central, East and Southeast Asian deformation. Their estimate agrees with earlier data (Huang et al., 2015). Based on Schellart et al.'s (2019) estimates, we have considered an indentation movement of 1,800-2,000 km in our model experiments. This India-Asia convergence requires an indentation rate, that is, the advance velocity of the leading edge of the Indian plate, ∼3.5 cm/yr, when averaged over ∼50 Ma (Schellart et al., 2019). The geological records, however, clearly suggest that the India-Asia collisional kinematics never remained steady, but decelerated with time (Copley et al., 2010;Molnar & Stock, 2009;Pusok & Stegman, 2020). We thus aimed to explore from our laboratory model experiments the possible effects of such decreasing indentation velocity on the crustal flow and deformation patterns in Tibet. To recognize this effect, we performed three types of experiments, and compared their results. Type 1 experiment simulated a steady fast-rate collisional kinematics, setting the indentation rate (I FAST ) at a constant value of 5.5 cm/yr. Type 2 experiment simulated a steady, but slow-rate collisional kinematics with a constant indentation rate (I SLOW ) of 3.5 cm/yr, which represents the average indentation rate. In the Type 3 experiment, we generated an unsteady collisional kinematics, beginning with a fast indentation rate of 5.5 cm/yr, but decelerating (I DEC ) to 3.5 cm/yr with time. This experiment exactly represents the reducing India-Asia convergence velocity reported in the literature (Copley et al., 2010;Molnar & Stock, 2009). In all the three types of experiments, the Indian indenter's leading edge was positioned at a distance around 1,800-2,000 km south from its present-day position. The total indentation of 1,800-2,000 km in the model is expressed in percentage as ϵ = 100%.

Type 1 Model Experiment: Steady Fast-Rate Collision
In this experiment of steady fast-rate (I FAST = 5.5. cm/yr) collision at ϵ = 30% convergence the model produced a NNE directed high velocity (5 cm/yr) crustal flow field in southern Tibet, gradually weakening nearly to 0 cm/yr in the northern margin of Tibet (Figure 3a-i). When ϵ = 50%, the NNE directed flow field grew northward, setting strong flow velocities (∼4 cm/yr) in the central Tibetan region. But, the northernmost part of the plateau still remained nearly undisturbed (Figure 3a-ii). At ϵ = 70%, the northward crustal flow interacted with the rigid Tarim basin situated in the north, and started to deflect towards north-eastern Tibet. The crustal flow velocities also showed a strong spatial heterogeneity, ∼5 cm/yr dominantly in southern and central Tibet, whereas ∼3 cm/yr in northern Tibet (Figure 3a-iii). This velocity gradient clearly suggests NNE-SSW contractional deformations in a large part of Tibet. At ϵ = 100%, that is, after the completion of total ∼2,000 km Indian indentation to the present-day condition, the entire Tibetan plateau continued to evolve with dominant NE-SW contractional deformations, maintaining crustal flow velocities nearly 5 cm/ yr in southern Tibet, ∼3.5 cm/yr in central Tibet, and ∼2.5 cm/yr in northern and north-eastern Tibet . During the entire stage of collisional event, the model did not produce any clockwise rotation of the crustal flows around EHS, as observed in the present-day Tibet ( Figure 1a).
We plotted the principal strain axes on the second invariant of strain-rate tensor maps to study the spatio-temporal variations of crustal deformation in Type 1 model experiments. The plots show intense deformation (high strain rates ∼4.5 × 10 −8 yr −1 ) localization preferentially near the leading edge of the Indian indenter in the initial stage (ϵ = 30%) of India-Asia collision (Figure 3b-i). With progressive collision the high strain-rate zone propagated northward. During ϵ = 30%-70% convergence the northern and the north-eastern regions of Tibet had crustal shortening at strain rates of 6 to 2.0 × 10 −8 yr −1 (Figures 3b-ii, iii, iv). The principal strain axes plot suggests that the entire collisional event underwent contractional deformations broadly in the NNE-SSW direction.
The topographic analysis reveals that at 30% convergence, the plateau developed a northward topographic slope of 0.02°, suggesting a higher uplift rate in southern Tibet, compared to the northern region. At this stage the plateau gained virtually no topographic slope (0.01°) in the east-west direction (Figure 3c-i). For ϵ = 70%, the northward slope steepened further to 0.1°, but the E-W topographic slope remains almost steady at 0.02° (Figure 3c-ii). Increasing convergence resulted in a topographic modification, marked by the increase of the northward topographic slope from 0.1° to 0.6° when ϵ = 100% (Figure 3c-iii).

Type 2 Model Experiment: Steady Slow-Rate Collision
The Type 2 experiment was similarly run at a steady, but relatively slow rate (I SLOW = 3.5 cm/yr) to cover the total ∼2,000 km convergence. The model run at ϵ = 32% produced high crustal flow velocities (3.5 cm/ yr) in the southern Tibetan region, where the flow had broadly a NNE-ward trend. However, central and northern Tibet were left virtually undisturbed, showing little or no crustal flow (Figure 4a-i). At ϵ = 50% the high-velocity (∼3 cm/yr) crustal flow regime in the southern region propagated northward to set a strong crustal flow (∼2 cm/yr) in central Tibet (Figure 4a-ii). The northern part of Tibet also became kinematically active, but with significantly low velocities (∼1 cm/yr). After ϵ = 70% strong high-velocity crustal flows still dominated in southern (∼3 cm/yr) and central Tibet (∼2 cm/yr), and at the same time the flows in northern Tibet -strengthened by increasing their velocity to ∼1.4 cm/yr (Figure 4a-iii). The rigid Tarim basin deflected crustal-flows towards north-eastward, but at much lower magnitudes, as compared to that observed in the Type 1 experiment. At the completion of experimental run (i.e., ϵ = 100%) the model maintained the same kinematic state of crustal flows with slightly reduced magnitudes (∼3 cm/yr in southern Tibet and ∼1.8 cm/yr in northern Tibet; Figure 4a-iv). This Type-2 experiment also failed to produce the clockwise turning crustal-flow field around EHS of the present-day Tibet (Figure 1a).
The model run at ϵ = 30% produced a gentle northward slope (0.02°), but negligibly small (0.009°) eastward slope (Figure 4c-i). With increasing convergence the northward topographic slope of Tibet steepened from 0.02° to 0.03° when ϵ = 70%. Meanwhile, the eastward topographic slope had a small value, ∼0.01° To highlight the main findings from the two steady-state collision experiments (Type 1 and Type 2), irrespective of fast or slow convergence rates, the collision did not produce eastward crustal flows or clockwise rotation of the crustal flow field around EHS in model Tibet at any stage of the convergence event. None of them gave rise to the present-day E-W extensional strain in central Tibet, as revealed from GPS observations (discussed later in section). We, therefore, chose to run the unsteady collision experiments with decelerating convergence rates, as applicable for the India-Asia collision event (Copley et al., 2010;Molnar & Stock, 2009).

Type 3 Model Experiment: Decelerating Collision
To perform this unsteady collision experiment, we reduced the convergence velocity (I DEC ) stepwise in three stages, approximated to the decreasing India-Asia convergence velocity shown in Figure 1d. In stage S1 (ϵ = 0%-40%; I DEC = 5.5 cm/yr), the model produced a strongly heterogeneous NNE-directed crustal-flow regime at ϵ = 30%. The flow occurred at high velocities (∼5 cm/yr) in southern Tibet, but decreased to almost zero at the northern boundary (Figure 5a-i). In stage S2 (ϵ = 41%-70%; I DEC = 4.5 cm/yr), the flow regime propagated northward to induce crustal flow velocity ∼4.2 cm/yr in central plateau at ϵ = 50%. The model produced flows also in northern Tibet, although at much lower rates (∼1.5 cm/yr) (Figure 5a-ii). At this stage the rigid Tarim basin acted as a mechanical barrier, and reduced the effective width of the western Tibetan plateau to make it much narrower than its eastern flank. This west to east variation in the mechanical setting caused the crustal flows to take a north-east ward turn in the eastern parts of Tibet (Figure 5a-iii).
In stage S3 (ϵ = 71%-100%; I DEC = 3.5 cm/yr) of collision, the model showed a remarkable change in the crustal flow pattern (Figure 5a-iv). Southern Tibet underwent NNE-directed flows with an average rate of 3.5 cm/yr, which reduced to ∼2 cm/yr in north-eastern Tibet with an azimuthal change grossly to ENE. Second, the flow field developed a distinct region of east-directed crustal flows with magnitudes of 1.5-2 cm/ yr in the central and southern parts of the plateau (Figure 5a-iv), which was not observed in the steady collision experiments. Interestingly, the eastward crustal flow, after facing resistance from the mechanically strong Sichuan basin in the east, took a clockwise turn to channelize through the passage in between the EHS and Sichuan basin, and escaped southward, forming a fan-like front. This region accelerated the crustal flows, with magnitudes reaching ∼3 cm/yr (Figure 5a-iv). At the final stage (ϵ = 100%), the model produced a crustal-flow pattern strikingly similar to the present-day GPS velocity field ( Figure 1a; Liang et al., 2013;Zhang et al., 2004).
The decelerating collision kinematics led to characteristic topographic evolution in model Tibet. In stage S1 (ϵ = 30%) the southern plateau gained higher topographic elevation to form a gentle northward slope of 0.03°, and almost an equal slope of 0.04°in the E-W direction (Figure 5c

A Comparative Analysis of Type 1, 2 and 3 Experiments
All the three types of model experiments showed that the rigid Tarim basin at the northern margin of Tibet acted as a mechanical barrier to the northward crustal-flows induced by the Indian indenter. The rigid block forced the crustal flows to redirect their path in the northeast direction (Figures 3a, 4a, and 5a), and contributed to the development of relatively higher topography in western Tibet. The three types of experiments also suggest that the deformation fronts in Tibet progressively migrated northward after the initial phase of collision (Figures 3b, 4b, and 5b). However, the spatio-temporal variations of crustal flow patterns and topographic upliftment produced in the three models are markedly different from one another. Type 1 and Type 2 models maintained a constant northward to north-eastward crustal-flow azimuth, for the entire period of experimental run. They never set any dominant eastward flow at any stage of the convergence history. In contrast, the Type 3 model produced a distinct eastward flow field after 70% convergence when the indentation velocity (I DEC ) reduced to 3.5 cm/yr (Figure 5a).

Eastward Crustal-Flow Pattern
We will now compare the present-day (for ϵ = 100%) crustal-flow patterns produced in the three types of model experiments with the present-day GPS velocity field ( Figure 6). Type 1 model produces crustal-flows dominantly in the NNE direction, which continued to occur with high velocities in the extreme north-east corner, without any sign of clock-wise turn against the EHS (Figure 3a-iv). The flow pattern in eastern Tibet under this steady fast-rate collision condition markedly disagrees with the GPS velocity field (Figures 6a  and 1a). Similarly, the model topography (steep northward slope, 0.6° and 0.04° E-W slope; Figure 3c-iii) does not match with the present-day topographic slopes of the Tibetan plateau (Figure 1b). Type 2 model also shows a discernible mismatch of the crustal-flow pattern with the GPS velocity field (Figure 6b). The high-velocity crustal flow regime occurs only in southern Tibet, but without much transmission of contractional deformation in north-eastern Tibet (Figure 4a-iv).
The crustal flow pattern produced in the Type 3 model shows an excellent agreement with the present-day GPS velocity field in Tibet (Figure 6c). The model topographic slopes along E-W and N-S also show a good match with the present-day topographic slopes (Figures 5c-iii, and 7b). Moreover, a time-series analysis of topographic uplift from the Type 3 model indicates that the rigid Tarim basin in the north acted as a mechanical obstruction to reduce the horizontal crustal flow velocity, and facilitated western Tibet to grow vertically faster at an average ∼0.24 mm/yr than the eastern Tibet at ∼0.10 mm/yr in the initial 25-30 Ma of the India-Asia collision. Such differential uplift rates resulted in a first-order eastward topographic slope during the phase (∼50-19 ± 3 Ma) of high indentation velocity (≥4.5 cm/yr) (Figure 7c). The topographic slope underwent gravity collapse in the phase of reduced indentation velocity to 3.5 cm/yr, which in turn set the eastward crustal flow in central Tibet. The gravitational collapse and onset of the eastward crustal flow can be theoretically interpreted in terms of the Argand number Ar, which is the ratio between two stresses, one arising from crustal thickness contrasts and the other required to deform the material at the ambient strain rates (England & McKenzie, 1982). According to this theoretical formulation, the balance between the gravity and the resistive viscous forces determines the rate of gravity-driven crustal flows in an orogen. A decrease in the convergence velocity results in an increase in Ar, that is, a greater influence of the gravitational potential energy (GPE) in the orogen dynamics. In the higher Ar condition gravity forces arising from the crustal thickness contrast thus captured the dynamics of Tibetan tectonics, setting in an eastward flow in the plateau when the convergence rate significantly dropped in Stage-S3. On the other end, the rigid Sichuan basin on the eastern flank of Tibet resisted the rapid eastward crustal-flow originated from the central plateau, and eventually forced it to take a clockwise turn around the EHS and flow towards the south with a fan-like front in the Burmese plate region (Figures 5a ii-iii). Topographic analysis suggests the Type 3 model developed a steep (0.1°) eastward slope after 70% convergence, which subsequently reduced to 0.04° at the present-day (Figure 5c). This clearly indicates that the eastward flow was driven by topography-controlled gravity current (Copley & Mckenzie, 2007), which occurred in the decelerating kinematic condition of Indian plate indentation. Earlier laboratory models for India-Asia collision also suggested that the gravitational collapse of the Tibetan plateau governed the late stage extensional tectonic in the orogen (Bajolet et al., 2013(Bajolet et al., , 2015Fournier et al., 2004). Fournier et al. (2004) showed that a confined western and a free or weakly confined eastern boundary of Tibet directed the gravitational spreading of Tibet in the east direction. Our model adds that the combined effect of lateral crustal heterogeneity (Tarim and Sichuan basins) and a reducing India-Asia convergence velocity controlled the differential topographic upliftment from west to east, and its subsequent gravitational collapse in the late stage (last ∼18 Ma). This model also suggests that the rigid Sichuan basin in the eastern boundary of Tibet forced the eastward flow to take a turn in the south eastern direction, forming a fan like front in the Burmese platelet.

Contractional to Extensional Strain Fields
The three types of experiments had remarkable differences in the spatio-temporal variations of their strainrate fields. In Type 1 and Type 2 models the contractional zones progressively migrated from southern to Comparison of the present-day model topographic profiles, obtained from the Type 3 model experiment (I DEC = 5.5 to 3.5 cm/yr) with the DEM profiles of the present-day plateau. (c) Calculated 3D topography from the Type 3 model produced after 70% of total convergence. Note that we have observed different rates (average) of topographic uplift from west to east in the Tibetan plateau.
northern Tibet, albeit with different magnitudes. The contraction direction was dominantly NNE-SSW to NE-SW. Both the models failed to produce E-W extensional strain in central Tibet (Figures 3b and 4b). In contrast, the Type 3 model at the late stage (i.e., after 70% of total convergence) gave rise to strong ESE-WNW extensional strains in central Tibet (Figure 5b-iv), which is a characteristic feature of the present-day deformation in Tibet ( Figure 8a). Interestingly, this last stage of convergence had contraction deformations in north-eastern Tibet and the region between EHS and Sichuan basin. However, the E-W contraction in the region between EHS and Sichuan basin showed a transition to E-W extension down to the southern region of Burmese platelet. We have compared the strain-rate field of the Type 3 model experiment with the present-day geodetic strain rates (Figures 8a and 8b). Their remarkable match suggests the Type-3 model experiment provides the best approximation to the kinematic history of the Tibetan plateau.

Orogen-Parallel Shear Deformations
We have developed a shear strain-rate     xy map from the Type 3 model velocity field for the present-day Tibetan plateau. The map shows shear partitioning in laterally persistent orogen-parallel belts across model Tibet. Strong sinistral shear motion (−1.5 × 10 −8 yr −1 ) localizes in northern Tibet, whereas dextral shear motion of nearly equal magnitude (1.5 × 10 −8 yr −1 ) concentrates in southern Tibet (Figure 8c). The transition between sinistral and dextral shear occurs in central Tibet, where the magnitude of shear rate becomes low (Figure 8c). The comparative analysis of the model results obtained from the three types of experiments suggest that the Type 3 model convincingly accounts for the effects of both lateral crustal heterogeneity and decreasing India-Asia indentation velocity in setting the present-day eastward crustal flow and the heterogeneous strain rate patterns in Tibet. We thus discuss the implications of our experimental study in the light of Type 3 model results in the following discussions.

Tectonic Evolution of the Tibetan Plateau: Model Versus Nature
Our laboratory model results suggest that the high-velocity crustal-flows localized initially close to the leading edge of the Indian indenter (Figure 5a), and resulted in rapid topographic uplift in southern and central Tibet (Figures 5c i-ii) during the initial stage of fast India-Asia collision (5.5-4.5 cm/yr). In the course of collision the high-velocity front propagated north to north-eastward (Figures 5a and 5b). This kinematic change led to northward migration of the topographic uplift in Tibet. However, in the late stage of collision with slow rates (3.5 cm/yr) north-eastern and eastern Tibet became active to develop elevated topography. On the other hand, the geological records used to establish the topographic evolution of the Tibetan plateau are not straightforward, but pose several debated issues. A line of evidences has been put forward to show that the central and some parts of southern Tibet already attained their elevation of 3-4 km by 40 Ma to form plateau topography, referred to as proto-Tibetan plateau (Kapp et al., 2005;Murphy et al., 1997;Y. Li, Wang, et al., 2015;Wang et al., 2014). Several workers have proposed the collision between Lhasa and Qiangtang block during the Late Cretaceous-Paleogene time (Jolivet, 2015;Kapp et al., 2005;Murphy et al., 1997), prior to the India-Asia collision at 50 Ma as the responsible factor for the growth of this pre-existing plateau. In counter to this view, (Spicer et al., 2021), based on paleontological proxies for paleoaltimetry, have recently shown central Tibet as a low-elevation region between the Gangdese Mountain in the south and the Qiangtang terrain in the north, which was originally below the sea level during the Cretaceous-Paleogene time. From magnetostratigraphic and radiochronologic dating Fang et al. (2020) have provided a similar interpretation for central Tibet topography with a low elevation (<2.3 km) at 39.5 Ma, which gained a height of 3.5-4.5 km at around 26 Ma. On the other hand, Y. Li, Wang, et al. (2015) compiled the available thermochronological data to demonstrate that the central Tibetan plateau, in overall, continued to grow till ∼20 Ma. Their analysis also stressed upon the fact that the Tibetan plateau progressively shifted its uplift location from the central to the northern region. This northward migration of the plateau upliftment is consistent with our model observations. The spatio-temporal variations of model crustal flows gave rise to heterogeneous strain distributions in the entire history of the Tibetan plateau (Figure 5b). In the stage (∼50 to 19 ± 3 Ma) of high convergence rates (5.5-4.5 cm/yr) the NNE-SSW contraction determined the tectonic deformation patterns in the plateau (Figures 5b i-iii). When the convergence rate was reduced to 3.5 cm/ yr (the present-day average value) at 19 ± 3 Ma the contractional mode of tectonic deformation started replacing by a distinct ESE-WNW extensional tectonics in southern and central Tibet (Figures 5b-iv and 8b) due to a topography-controlled gravity current (Copley & Mckenzie, 2007;Liang et al., 2013). This tectonic transition predicted from our scaled laboratory experiments is consistent with the deformation history of the Tibetan plateau. Geological imprints suggest that the E-W striking thrust faults in central Tibet were mostly active in the period of ∼50 to ∼22 Ma, and they subsequently became dormant. Active thrust faulting then shifted to localize preferentially in the southern and north-eastern margins (Y. Li, Wang, et al., 2015;Wang et al., 2014). After 19 ± 3 Ma the onset of ESE-WNW extensional tectonics manifested in the form of N-S striking, low-angle normal faulting in central and southern Tibet (Ge et al., 2015;Kali et al., 2010;Lee et al., 2011;Y. Li, Wang, et al., 2015;McCallister et al., 2014;Styron et al., 2013). Most of the available field geological data suggest that this transition took place across the orogen over a timespan, broadly in middle Miocene. For example, the N-S stretching extensional grabens and normal faults were initiated in southern Tibet at 14-13 Ma (McCallister et al., 2014), ∼13-12 Ma (Kali et al., 2010;Lee et al., 2011), ∼11 Ma (Garzione et al., 2003), <12.5 Ma (Edwards & Harrison, 1997) etc. In further north, similar ages are reported, 16-12 Ma (Styron et al., 2013), 11-8 Ma (Woodruff et al., 2013), ∼8 Ma (Harrison et al., 1995) and in northern Tibet it is ∼13.5 Ma (Blisniuk et al., 2001). Recently,  have suggested monotonically decreasing ages of rift developments across the orogen (from ca. 23 Ma at Leo Pargil to 2-3 Ma in the east at Cona). This observation is consistent with our present tectonic model, which produced GPE gradient-driven eastward lithospheric flows with reducing convergence rates. However, some of the recent studies have reported present-day thrust activities from Gangdese (Taylor et al., 2021) and Qimen Tagh areas Yakovlev, 2015) of southern and central Tibet. These observations suggest that the thrust systems in central Tibet, although they are not major contributors to the ongoing Tibetan tectonics, are not totally inactive.
Our experimental model in the late stage, that is, from 19 ± 3 Ma to the present day, produced a strong E-W crustal flows (∼1.5 cm/yr) in central Tibet (Figure 8b). The model estimate suggests an extensional strain of 0.8 to 2.2 × 10 −8 yr −1 in this region (Figure 8b-ii), which fairly agrees with the geodetic estimates (1.6 to 2.2 × 10 −8 yr −1 , Ge et al., 2015). The time series analyses of our model crustal-flow and the strain pattern support a two-stage evolution of the Tibetan plateau: (a) a prolonged period of contraction since the beginning of collision to as late as 19 ± 3 Ma, which was coupled with (b) a phase of wide-spread horizontal extension in the mid-Miocene time (19 ± 3 Ma) preferentially in central Tibet Maiti & Mandal, 2021;Wang et al., 2014).

Contractional Deformations in the NE Flank of Tibet
The model experiments suggest that north-eastern and eastern Tibet in adjacent to the Sichuan basins localized southward crustal flows at high rates (∼3.5 cm/yr) in the late stage (i.e., after 70% convergence) of India-Asia collision after ∼22 Ma (Figures 5a iii-iv). Model topographic analysis reveals that these regions gained significant topographic elevations after 19 ± 3 Ma (Figure 5c-iii). The strain-rate plot also indicates NE-SW contractional deformations in north-eastern Tibet and an E-W contraction near Sichuan basins during this late stage (Figures 5b iii-iv). Using thermochronological data several workers have constrained the initiation of upliftment/exhumation of the Longmen Shan thrust (LMST, Figure 7a, boundary between Sichuan basin and eastern Tibet) at around 22 Ma (Wilson & Fowler, 2011;Xu & Kamp, 2000). Although some studies indicated the existence of early elevated topography in Longmen Shan (X. Chen et al., 2003;Wallis et al., 2003), it has been suggested from local structural analyses and river incision data that the rapid exhumation and high topographic elevation occurred in this region dominantly during 12-10 Ma (Clark et al., 2005;E. Wang et al., 2012;Xu & Kamp, 2000). Moreover, some thrusts (i.e., Ganzi-Litang thrust belt, GLT) in the western part of the Longmen Shan region showed rapid exhumation during 20-16 Ma (Figure 7a). The deformation history of these thrust belts show an eastward migration of their activities from GLT in the west to Longmen Shan in the east. All these data reveal a Miocene uplift event in the eastern-most Tibetan region, which strongly supports our model results. On the other hand, some studies suggested the Qaidam basin in north-eastern Tibet to have developed its elevated topography during 40 to 22 Ma E. Wang et al., 2012;Yin et al., 2008). However, the existence of such elevated early topography is still debated (Kent-Corson et al., 2009;L. Wang, Cheng, et al., 2021). On further north-east, in the Qilian Shan region several NE-directed thrust faults (i.e., Nanshan Thrust belt, NST; Eastern Qlianshan-Nanshan thrust belt, EQNT; Northern Qilian Shan thrust belt, NQLT) became active in between 33 and 9 Ma (Yin et al., 2002). Among them, NQLT is still active at present day (S. Yang et al., 2007; Figure 7a). These observations suggest that the major uplift in the NE region occurred in between 24 and 9 Ma, although some parts of Qilian Shan region started to gain their elevation since 40 Ma Bian, Gong, Chen, et al., 2020;Clark et al., 2010;B. Li et al., 2020;L. Wang, Cheng, et al., 2021). Our model also indicates that a small amount of the contractional deformation propagated northward to affect north-eastern Tibet during the fast-collision period (Figures 5a-ii, 5b-ii), but the major contractional deformations in the extreme north-eastern margin of Tibet was a late (19 ± 3 Ma to present-day) tectonic event (Figures 5b iii-iv). Several studies, however, have shown a complex pattern of topographic uplift in Tibet due to the presence of numerous pre-collision arcs and mechanically weaker lithospheric segments, as well as reactivation of suture zones (Zuza & Yin, 2017). For example, numerical modeling of Bian, Gong, Chen, et al. (2020) suggests that Qilian Shan region underwent early Cenozoic uplift preferentially as a consequence of its weaker lithospheric rheology. L. Wang, Cheng, et al. (2021) have recently estimated the early Cenozoic topographic elevation of Qilian Shan varying between 0.4 and 1.5 km, which is relatively higher than that in Eastern Kunlun Shan (0.4-1 km) located to its south. However, during Neogene time Eastern Kunlun Shan had tectonic uplift more than Qilian Shan (located in far north-eastern Tibet). With the present-day topographic elevation of ∼3-4 km in Eastern Kunlun Shan and ∼2-3 km in Qilian Shan, it is clear that, despite the occurrence of some pre-existing topographic relief in north-eastern Tibet the topographic upliftment progressed from south to north, where GPE controlled north-east wards crustal flows might have played an important role, as inferred by L. Wang, Cheng, et al. (2021).
The present-day NE-SW contractions in north-eastern Tibet occurs at a strain-rate of 2 to 3.2 × 10 −8 yr −1 , estimated from GPS data (Gan et al., 2007). Our model estimate yields a similar order of strain-rates (2 × 10 −8 yr −1 ). The geodetic observation also indicates that the contraction rates in the region between EHS and Sichuan basin are measured as 3 to 4 × 10 −8 yr −1 (Gan et al., 2007), which fairly matches with our model estimates, 3.2 to 4.2 × 10 −8 yr −1 (Figures 8a and 8b). Based on the correlation between the model and natural data, we suggest that the contractional deformation in NE Tibet is driven by the NE-directed crustal flows from the plateau interior. On the other hand, in eastern Tibet the eastward crustal-flow encountered a strong resistance from the Sichuan basin to localize a contractional strain field. Clark and Royden (2000) proposed a mechanical model to explain the development of topography in eastern Tibet. According to their model, the topographic pressure driven lower crustal flows determined the topographic elevation. Our model also suggests GPE controlled flows, but on a lithospheric scale, as the underlying factors for topographic growth in eastern and north-eastern Tibet.

The Present-Day Crustal Velocity Field and the Plateau Topography
Type-3 model estimates yield the present-day velocity magnitudes and their directions in southern (∼3.5 cm/ yr), central (∼3.2 cm/yr), and north-eastern Tibet (0.5-1 cm/yr) (Figures 5a-iv, 6c), which match well with the GPS velocities, ∼3.8 cm/yr, ∼3.5 cm/yr, and 0.5-1.5 cm/yr in the respective regions (Gan et al., 2007;Liang et al., 2013;Zhang et al., 2004). In addition, the model reproduces a spectacular clockwise rotation of the velocity vectors (1-2 cm/yr) around EHS, as predicted from GPS studies (Zhang et al., 2004). The present-day topographic gradients in the north-south (AB), east-west (CD) sections as well as in north-eastern and south-eastern model Tibet also show a good match with those obtained from the digital elevation maps of Tibet (Figure 7b). The finding validates our model interpretation that the eastward first-order topographic gradient drives the present-day west to east crustal flows, as shown in earlier models (Cook & Royden, 2008;Copley & Mckenzie, 2007;Klemperer, 2006;Liu & Yang, 2003).

Eastward "Glacier-Like" Flow and Its Implications
We calculated the E-W component of crustal-flow velocities from the Type 3 model that represents the present-day tectonics of Tibet. Its N-S variation describes a parabolic velocity profile with the maximum magnitude located in the central part and the minimum in the northern and southern margins of Tibet ( Figure 9a). The parabolic velocity profile typically represents a Poiseuille type flow, implying that the eastward flow is driven by a pressure gradient (Yin & Taylor, 2011). Yin and Taylor (2011) showed that, in the last 15 Ma the continental lithosphere (on a horizontal scale of 1,000 km) in central Tibet flowed eastward like a fluid mass, between two rigid walls, relatively moving towards each other (Bajolet et al., 2015;Fournier et al., 2004). They compared their model observations with existing GPS data, and proposed horizontal shear applied at the base of the Tibetan mantle lithosphere or lateral spreading of the thickened Tibetan lithosphere as a driving factor for the late Cenozoic deformation in central Tibet. Our model suggests that the first-order E-W topographic gradient induced a pressure difference, which eventually acted as the driving force for this flow, as proposed by Copley and Mckenzie (2007) for the gravitational spreading of crustal materials in southern and south-eastern Tibet. The eastward-flow began only after the indentation velocity decreased to the present-day average of 3.5 cm/yr (Figures 5a-iii-iv). A comparison of the model E-W velocity fields with the GPS estimates of Zhang et al. (2004), also suggest that their parabolic profile is a characteristic representation of the present-day kinematics in the central Tibetan plateau (Figure 9b). The model experiments reveal another interesting fact that the eastward flows strengthen their velocity magnitudes in the east direction, 0.5 cm/yr to nearly 2 cm/yr, which conforms to GPS observations (Figure 9c). Both the model and GPS data show the north-directed crustal flows with decreasing magnitudes from south to north, 3.5 cm/yr to ∼1 cm/yr (Figure 9d).
Our model results discussed above suggest that the decelerating collision kinematics set a Poiseuille-type eastward crustal-flow in central Tibet, where the Tarim basin in the north and the Himalaya in the south acted as the guiding walls (Figures 9a and 9b). This horizontal flow resulted in sinistral shear strain in the north and dextral shear strain in southern Tibet (Figure 8c). The two shear fields of opposing shear senses explain most of the strike-slip faults (Taylor & Yin, 2009) and their motion in southern and northern Tibet (Figure 8c). Yin and Taylor (2011) combined field observations with their mechanical models (localized basal shear or lateral gravity spreading) to demonstrate the kinematics of these two shear domains, and explained the occurrence of right and left strike-slip faults in central Tibet as a conjugate strike-slip fault system (Figure 8c). Their model, however, does not account for the origin of non-conjugate right-and left-slip faults on either side of the conjugate fault system. Our model Tibet, on the other hand, produces two shear fields in the extreme southern and northern parts of Tibet, with a transition zone of little or no shear in the central part. We suggest that the conjugate faults with opposing slip senses have formed in this transition region dominated by symmetrical eastward flows, whereas non-conjugate left-and right-slip strike faults in the shear dominated fields of its northern and southern sides, respectively.

Role of Lateral Crustal Heterogeneities in Tibetan Tectonics
In the laboratory model experiments the rigid Tarim basin lying on north of Tibet acted as a mechanical barrier to the northward crustal flows in the Tibetan plateau, and redirected the flows in the NE direction (Figures 10a and 10b). Moreover, this resistance forced the Tibetan crust to uplift at higher rates in its western  Zhang et al. (2004) along the green, red and blue lines (marked in Figure 1a). (c) Comparison of the trend of increasing eastward velocity between experiment and GPS data along CD (Figure 1c). (d) Decreasing northward velocity in present-day Tibet, obtained from experiment and GPS data along the green, red and blue lines (marked in Figure 1a). side, relative to its eastern counterpart (Figure 7c). Such a differential uplift eventually developed first-order topography with an eastward surface slope. Interestingly, Tibet still maintains an eastward slope of 0.03° (Figure 7b). Numerical models of Yang and Liu (2013) also produced results similar to that observed in our viscous analogue model, where Tarim basin in the north produced higher topographic upliftment in the western Tibet, compared to east. Their model also predicted that the zone of crustal uplift propagated from south to north in the early stages and west to east in the later stages. However, their model Tibetan plateau attained a maximum topographic elevation of 5 km in the last 5-0 Ma of India-Asia collision. In contrast, the same maximum elevation was attained in our model by 19 ± 3 Ma, as also estimated from various field observations (DeCelles et al., 2007;Fang et al., 2020).
Our model experiment strongly supports the hypothesis proposed by several earlier workers that the eastward flows encountered a resistance from the mechanically strong Sichuan basin in eastern Tibet (Clark et al., 2005;Cook & Royden, 2008;Penney & Copley, 2021). The Sichuan basin thus acted as a crucial factor in forcing the crustal flows to take a clockwise turn around the EHS (Figure 10c). Interestingly, this experimental finding holds an excellent agreement with previous numerical India-Asia collision models (Cook & Royden, 2008;Pusok & Kaus, 2015;Yang & Liu, 2013). In addition, the present model allows us to propose that the fan-like flow divergence in the Burmese plate region occurred dominantly in a later stage (after ∼19 Ma) of the India-Asia collision when the convergence rate reduced to 3.5 cm/yr. encountered resistance from rigid Tarim basin, and were deflected towards northeast. At this stage the western and central Tibet uplifted at relatively higher rates, compared to eastern Tibet, creating an eastward topographic gradient (0.1°). (c) At 19 ± 3 Ma (ϵ = 70%), decrease of I DEC to the present average (3.5 cm/ yr), resulted in gravitational collapse of the plateau. The eastward topographic gradient directed the crustal flow towards east to initiate the E-W extensional tectonics, as observed in the present-day configuration. Note that, during this stage the eastward crustal flow followed a Poiseuille type velocity profile, giving rise to sinistral and dextral shear in northern and southern Tibet, respectively.
The interaction of Tarim basin with the Tibetan plateau in the course of India-Asia collision has been quite complex, showing both underthrusting of Tarim and strike-slip faulting (S. Wang, Chen, et al., 2021;Wei et al., 2013). From teleseismic studies Kao et al. (2001) investigated lateral structural changes around Tarim and inferred that the basin behaved to be rigid with little shortening. They also found the Moho discontinuity of Tarim dips gently toward the south to ∼50 km under the Kunlun foreland, and proposed for a complex pattern of lithospheric collision than the simple configuration of one block subducting beneath another. By correlating seismic reflectors with the corresponding dated outcropping sections, S. Wang, Chen, et al. (2021) have shown that rigid Tarim lithosphere is underthrusting beneath the western Kunlun Shan region of Tibet since ∼30 Ma at an average rate of ∼8.18 mm/yr. This rate comprises shortening rate of ∼3.8-5.7 mm/yr and latitudinal propagation (northward) rate of the western Kunlun Shan ∼4.4-2.5 mm/ yr. In our experiments the model Tibetan crustal masses had a tendency to flow over the Tarim block at the late stage (after 19 ± 3 Ma), suggesting the Tarim block relatively underthrusting against the Tibet ( Figure  S4-b in Supporting Information S1).

Limitations and Future Scopes
Our laboratory modeling considers a two-layer mechanical structure of the Himalaya-Tibet system, consisting of a top high-viscosity layer (crustal lithosphere) upon a relatively low-viscosity layer (mantle region). The model crustal lithosphere is treated as a single viscous sheet, scaled to represent the depth-averaged density and viscosity of a 40 km thick crust. The model thus reproduces the overall continental-scale deformation in Tibet and its first-order topography produced during collision, as in the case of earlier TVS approximations (England & McKenzie, 1982). This model excludes higher-order discontinuities (i.e., faults, folds etc.) and their topographic effects, which is a limitation of the present study. In addition, we have entirely constrained the model boundary conditions kinematically to incorporate time dependent variations in India-Asia convergence velocities based on the available geological records. However, changes in plate motion are necessarily driven by changes in one or more driving or resisting forces. Iaffaldano et al. (2006) showed that, for large mountain belts and plateaus the topographic loads consume a significant amount of the driving tectonic forces by increasing frictional resistance between the overriding and subducting plates. The changing force-balance condition can thereby govern the long-term evolution of plate motion. Our present model study, however, does not account for the dynamics responsible for reducing convergence rate in the course of the India-Asia collision, which is also a limitation.
Some workers suggested that rigid Tarim block at the northern margin and the Sichuan basin in the eastern flank are pushed northward and eastward, respectively in the course of the India-Asia collision (Pan & Shen, 2017;S. Wang, Chen, et al., 2021). In fact, the GPS observations support this hypothesis (Liang et al., 2013;Zhang et al., 2004). However, there is a considerable uncertainty with the exact timing of such movement. Moreover, Meyer et al. (1998) showed that the northern boundary of Tibet adjacent to the Altyn Tagh Fault region remained stable since the beginning of the India-Asia collision. Clark et al. (2010) have also shown that the northern and north-eastern margins of Tibet remained nearly at the same position relative to stable Eurasia in the entire period of India-Asia collision, barring a little northward movement. Considering this kinematics of the northern margin on geological time scales, we have chosen to impose a fixed-wall condition for these rigid blocks. There is a scope of the present experimental study to explore the effects of their present-day movement, as predicted from GPS velocities, on the crustal flow pattern in Tibet.
Our experimental models reproduce the horizontal strain rate tensors grossly in agreement with the geodetic strain rates in Tibet. However, several studies have shown contraction deformations at the southern margin of Tibet as a manifestation of transient elastic strain accumulation around the Himalayan thrust faults, that is, thrusting on the range margin (e.g., Stevens & Avouac, 2015). These transient elastic effects are not accounted for in our experiments, which is evidently a limitation in our study. The present model thus does not fully capture the governing dynamics of Tibet and some deformational observations, such as dominantly ESE-WNW permanent extensional strains with little or no NE-SW contraction in southern Tibet, as revealed from earthquake focal mechanism solutions (Armijo et al., 1986;Copley et al., 2011;Elliott et al., 2010;Styron et al., 2015). Also, our model Tibet considers a flat initial topography. However, there is an ample evidence for the existence of pre-collisional topographic elevation in central (i.e., Lhasa and Qiantang blocks) and north-eastern Tibet (Qilian Shan), as discussed earlier. These high-topographic zones are expected to undergo uplift at lower rates due to their greater GPE during the collision. Consequently, crustal shortening would preferentially localize more intensely in the low-lying areas. There is a scope for future studies to validate this hypothesis experimentally by incorporating pre-existing topography and crustal heterogeneities in the initial model setup. Despite all these limitations, our Tibet model successfully demonstrates the development of the first-order E-W Tibetan topographic slopes, and their subsequent collapse with decreasing Indian convergence rates, setting in eastward crustal flow in the late phase of India-Asia collision (Figures 7b and 10). An agreement of the model predicted present-day crustal-flow velocities ( Figure 6c) and the overall strain tensor fields (Figures 8a and 8b) with actual natural data provide additional strengths to this study.
In a recent study (Webb et al., 2017) it has been shown that the change in slab dynamics, anchoring of the Indian continental slab during 30 to 25 Ma and subsequent slab breakoff in the last 13 Ma, has critically influenced the Himalayan tectonics. The slab breakoff event resulted in a shift of crustal shortening from deep inside the orogenic wedge in the early phase to a rapid rise of the Mountain in the last 13 Ma years. In the present model we have considered a simple steady underthrust geometry, without any slab steepening, break-off, or flattening in the experimental setup. Recent numerical models of Dasgupta et al. (2021) allow us to stipulate that unsteady slab-dip (i.e., steepening or flattening) would be an additional factor in determining the contractional versus extensional tectonics of the overriding plate. Thus, this issue opens up a further scope for advancing the present Tibet model with an unsteady slab underthrusting.

Conclusions
The laboratory model experiments lead us to the following conclusions. (a) Steady-state kinematics of India-Asia collision at either high or low convergence rates would not produce the present-day crustal flow and deformation patterns in the Tibetan plateau, as predicted from GPS observations. (b) We interpret the onset of eastward crustal flow and related ESE-WNW extension in Tibet as a consequence of decelerating Indian indentation velocity (I DEC ), from ∼5.5 cm/year at ∼50 Ma to ∼3.5 cm/year at 19 ± 3 Ma in the India-Asia collision history. (c) The rigid Tarim block, located in the northern boundary of Tibet has played a critical role (Figure 10) in developing an eastward first-order topographic gradient by 19 ± 3 Ma (Figure 7c). This topographic gradient later dictated the crustal material to flow in the east direction, resulting in a transition of NNE-SSW contraction to ESE-WNW extensional tectonics and formation of widespread nearly N-S trending grabens in southern and central Tibet when I DEC decreased to ∼3.5 cm/yr at 19 ± 3 Ma. (d) The present model also suggests that the eastward crustal flow eventually encountered the rigid Sichuan basin in eastern Tibet (Figure 10c), which forced the flow to take a turn (clockwise rotation) in the southern direction. (e) The orogen-parallel velocity component in our model Tibet follows approximately a parabolic variation from its southern to northern boundary, which broadly agrees with the GPS velocity profile. We thus propose Poiseuille-type kinematics for the orogen-parallel flow component, driven by a pressure difference due to the west to east topographic gradient. This Poiseuille-type flow explains the localization of sinistral and dextral strike-slip faults preferentially in the northern and the southern parts of Tibet. Finally, our laboratory model estimates of extensional and contractional strain rates: 0.8 to 2.2 × 10 −8 yr −1 and 3.2 to 4.2 × 10 −8 yr −1 confirm the GPS-based strain rates in central Tibet and its eastern flank, respectively.