Effect of rheological heterogeneities on the lithospheric deformation of the Tibetan Plateau and neighbouring regions

The Tibetan Plateau, induced by the India-Eurasian collision, has the highest average altitude in the world. During its uplift vertically, the Tibetan Plateau has been considered to expand laterally. However, there are several strong and almost non-deformable cratons on its periphery, such as the Tarim, North China craton, and South China block. The present landform features show that these cratons limit the expansion of the Tibetan Plateau. However, there is still much controversy over whether the deformation can be transmitted to periphery orogens or reactivate ancient orogens in the cratons. This study used numerical models to investigate the effect of rheological heterogeneities on the lithospheric deformation of the Tibetan Plateau and its neighbouring regions. The results show that the lateral heterogeneities of the lithosphere have an important influence on the deformation or strain partitioning. Generally, during the lateral expansion of the Tibetan Plateau, its peripheral cratons can transmit the deformation or high strain to neighbouring weak orogens. This case can be used to understand the Tian Shan orogen, which was reactivated by the India-Eurasian collision. However, when the orogens inside the cratons have high lithospheric strength, high strain is difficult to distribute on them and the expanding Tibetan Plateau is constrained by its peripheral cratons. These results can be used to explain the ancient orogens that are not strongly deformed, such as the Jiangnan orogen in the South China block. Because these orogens formed at the same time as the cratons and have relatively high lithospheric strength. In addition, the large lithospheric thickness difference and low crustal rheological contrast favor high strain rates localized on the lithosphere of the ancient orogen in the craton, such as the Trans-North China orogen in the North China craton.

involved in the plateau (Tapponnier et al., 2001;Wang et al., 2014a;Yin and Harrison, 2000). The Tarim, North China craton, and South China block are the three main cratons in East Asia. However, these three cratons were not integrated into the plateau, on the contrary, they severely restricted the expansion of the Tibetan Plateau based on the landform features ( Figure 1A). There are three main orogenic belts outside or within these cratons: the Tian Shan orogen outside the Tarim, the Jiangnan orogen within the South China block, and the Trans-North China orogen within the North China craton. Whether the deformation or strain caused by the India-Eurasian collision can be transferred by cratons to these orogens is still controversial (Molnar and Tapponnier, 1977;Sun and Liu, 2018;Yuan et al., 2013).
Previous numerical studies have shown that during the lateral expansion of the Tibetan Plateau, if a stable craton is encountered, strain can be transferred to the weak orogen on the other side through the craton (England and Houseman, 1985;Neil and Houseman, 1997). For example, the Tibetan Plateau encountered the Tarim during the expansion of its northern margin. The strong Tarim, while undergoing little internal deformation, transfers strain to the Tian Shan orogen, producing significant crustal thickening ( Figure 1B). The Tian Shan orogen belongs to the Central Asian orogenic belts, which were formed in the Paleozoic (Windley et al., 2007;Xiao et al., 2013). The expansion of the Tibetan Plateau or India-Eurasian collision reactivated the Tian Shan orogen (Jolivet et al., 2010;Li et al., 2022;Tapponnier and Molnar, 1979). Similar tectonics exist in the eastern and northeastern expansion of the Tibetan Plateau. Across the eastern margin of the Tibetan Plateau and the Sichuan basin, there is the Jiangnan orogen within the South China block ( Figure 1C). It was formed in the Neoproterozoic due to the collision of the Yangtze craton and the Cathaysia block (Li et al., 1995). Across the northeastern margin of the Tibetan Plateau and the Ordos, the Trans-North China orogen is located within the North China craton ( Figure 1D). It was formed in the Paleoproterozoic due to the collision of the west (Alashan and Ordos blocks) and east (North China plain) blocks of the North China craton (Zhao et al., 2000;Zhao et al., 2005). Although these two orogenic belts are located inside the cratons, their lithosphere has undergone a strong post-modification since the Cambrian. Especially, in the Mesozoic, the lithospheric thinning occurred in the east of North China craton and South China block (Zhu and Xu, 2019). Both for the South China block and the North China craton, there were significant differences in the thickness of the lithosphere between the east and west parts. In particular, the thickness of the mantle lithosphere in the eastern North China craton has decreased significantly, from~200 km tõ 80 km (Chen et al., 2009;Zhang et al., 2019;Zhu et al., 2011), and large-scale magmatic activities have occurred. However, the Ordos in the west of the North China craton and the Sichuan basin in the west of the South China block both maintained the stability of the craton  (An and Shi, 2006). JNO: Jiangnan orogen; TNCO: Trans-North China orogen; NCP: North China plain.

Frontiers in Earth Science
frontiersin.org 02 without significant deformation. In general, older orogens have stronger or more stable lithospheres due to the cooling of the lithosphere (Hu et al., 2000). However, under the action of strong neotectonic stress, some intraplate rheological weak zones may be reactivated, such as the Sprigg's orogen in southeastern Australia (Dyksterhuis and Müller, 2008) and the Laramide orogen in the United States (English and Johnston, 2004). Compared to the Tian Shan orogen, it is not clear whether the expansion of the Tibetan Plateau can affect these two ancient orogenic belts through the strong Sichuan basin and Ordos, respectively (Zhao and Zheng, 2007).
Previous studies have shown that the far-field effect of the India-Eurasian collision has extended to the Altai Shan and east China based on the active tectonics and topography (Molnar and Tapponnier, 1975;Molnar and Tapponnier, 1977). During this process, high topography and strain are mainly distributed along the preexisting faults, which implied the first-order importance of the weaknesses (Kong et al., 1997). However, these studies focused on crustal stress transmission (Dayem et al., 2009;England and Houseman, 1986;Kong et al., 1997). The dynamic effects of the India-Eurasian collision may have extended much deeper. On the one hand, recent study shows the subducting Indian mantle lithosphere beneath the western Tibetan Plateau collided with the Tarim lithosphere controls the uplift of the Tian Shan orogen (Huangfu et al., 2021). This implies the importance of the rheology of the mantle lithosphere. On the other hand, tomography studies have shown that continuous low-velocity asthenospheric mantle structure extends from the Tibetan Plateau to east China at the depth of 200-300 km (Huang et al., 2009;Liu et al., 2004). Numerical simulations suggest that significant lateral extrusion of the asthenospheric mantle driven by India-Eurasian collision is possible in the presence of a low-viscosity asthenosphere (Liu et al., 2004). Therefore, one possible explanation for the Cenozoic rifting and volcanism in east China is mantle flow driven by India-Eurasian collision. These studies show that the upper mantle dynamics should be considered in investigating the deformation distribution during the India-Eurasian collision.
In this study, a two-dimensional thermo-mechanical finite element model was used to consider the lithospheric strength heterogeneities of the Tibetan Plateau and neighbouring regions. We tried to explore how deformation or strain distribution when the expanding Tibetan Plateau collides with peripheral cratons. Especially, by considering the lithospheric deformation of the craton and asthenospheric flow, we try to investigate if the ancient orogenic belts in the craton, such as the Jiangnan orogen or the Trans-North China orogen, can be reactivated like the Tian Shan orogen. The results would help us to understand the far-field effects of the India-Eurasian collision and the interaction of different dynamic factors.

Governing equations
Two-dimensional thermo-mechanical finite element codes have been developed based on previous studies (Sun et al., 2019;Liu, 2018, 2023). For the mechanical deformation, the visco-plastic rheology was adopted in the numerical model, and viscosity depended on the temperature and strain rate. The Arbitrary Lagrange-Euler (Fullsack, 1995) and the marker-in-cell method (Gerya and Yuen, 2003) were used to solve the coupling equations.
The following equations are used in the two-dimensional thermal-mechanical coupled finite element model. The momentum conservation equation is where P − σij 3 is pressure and σ ij is stress tensor. δ ij is the Kronecker delta tensor. σ ij −Pδ ij + τ ij and τ ij is deviatoric stress tensor. x i is spatial coordinates (in the 2D case employed here: x 1 x, x 2 y). ρ is density, g j [0, −g] T is the gravitational acceleration. Here we use the index notation, which implies summation over its range when an index is repeated in a term. ρ ρ 0 [1 − α(T − T 0 )] is density, depending on temperature T. ρ 0 is the reference density under the condition of T 0 298 K; α is the thermal expansion parameter. α 3 × 10 −5 K −1 is used in this study.
The mass conservation or incompressible continuity equation is where v j is velocity vector. The constitutive equation is where η ef f is effective viscosity. The strain rate tensor _ ε ij can be written as: The ductile behavior is driven from experimental uniaxial creep law: A is material constant, n is power law exponent, E is the activation energy, R 8.314 J/(mol · K) is Boltzman's gas constant and T is the temperature in K. Frictional-plastic deformation would occur if the stress state satisfied a yielding criterion. Here we ignore the effect of pressure and active volume, which are of minor importance for the upper mantle (Ranalli, 1997a). In this study, we used a pressure-depended Drucker-Prager yield criterion: where σ II (J 2 ) 1 2 ( 1 2 τ ij τ ij ) 1 2 is the second invariant of the deviatoric stress tensor, C 0 is the cohension, α is the angle of internal friction. The plastic viscosity for plastic flow is in the regions that satisfy the yield criterion. Here To quantify the intensity and state of deformation (extension or compression), we define the maximum shear strain rate γ and dilatation strain rate θ to study the deformation.
where _ ε max and _ ε min are the maximum and the minimum principal strain rates, respectively. γ indicates the intensity of deformation. Positive θ indicates that the deformation is dominated by elongation or extension, while negative θ implies compressive deformation.
The energy conservation equation is given based on the extended Boussinesq approximation (EBA) (Christensen and Yuen, 1985;Hansen and Yuen, 2000) ρC P zT zt where ρ is density, C P is heat capacity, T is temperature, k is thermal conductivity, v j is velocity vector, Q r is radioactive heat production, Q s τ ij _ ε ij is frictional shear heating, and Q a is the adiabatic heating term. The adiabatic heating term can be simplified by neglecting deviations of the dynamic pressure gradients. That is Q a αT DP Dt ≈ − αTρgv y .

Finite element mesh and boundary conditions
The constructed numerical model is 2,200 km in length and 400 km in depth ( Figure 2). From left to right, the initial width of the Tibetan Plateau and craton in the numerical model are 1,050 km and 1,150 km, respectively. In the vertical direction, the model is divided into 5 layers: upper crust, middle crust, lower crust, mantle lithosphere, and asthenosphere. The upper, middle, and lower crustal thicknesses of the Tibetan Plateau were set as 15.8 km, 18 km, and 18 km, respectively. For the craton, the upper, middle, and lower crustal thicknesses were set as 10 km, 15 km, and 15 km, respectively. The initial topography of the Tibetan Plateau and craton are 1.8 km and 0, respectively. The initial crustal structure satisfies isostatic equilibrium. For the reference model (case R), the initial lithospheric thicknesses of the Tibetan Plateau and craton are 100 km and 130 km, respectively. The initial thicknesses of different layers change linearly from the Tibetan Plateau to the craton.
The mechanical boundary conditions of the reference model (case R) are imposed as follows. A 20 mm/yr constant velocity is imposed on the left side of the model. The bottom and left sides of the model are freeslip boundary conditions. The top side of the model is a full free boundary condition. In this case, the top surface may move vertically because of mass conservation (or incompressibility). We set the topography at the center of the craton as the reference point (0 km), then corrected the y-axis coordinates of all markers, which is relative to the elevation of the craton ( Figure 2). For this reason, the topography shown in this study indicates a relative change, not an absolute elevation. The thermal boundary conditions are imposed as follows. The thermal boundary conditions are zT zx 0 on the left and right sides. The temperature imposed on the top and bottom of the model is 0°C

FIGURE 2
The two-dimensional finite element model and boundary conditions (A) and the initial effective viscosity of the reference model (case R) (B). Different colors in (A) indicate different materials. Units 1-4 are the upper crust, middle crust, lower crust, and mantle lithosphere of the Tibetan Plateau or a weak orogen in other cases, respectively; units 5-8 are the upper crust, middle crust, lower crust, and mantle lithosphere of the craton or a strong orogen in the craton, respectively (Table 1); unit 9 is the asthenosphere. White lines in (A) indicate the initial isotherm.

Frontiers in Earth Science
frontiersin.org 04 and 1,500°C, respectively. The initial temperature from the top to the bottom of the lithosphere increased linearly from 0°C to 1,350°C. From the bottom of the lithosphere to the bottom of the model, the temperature increases linearly from 1,350°C to 1,500°C (400 km depth).

Parameters
Cenozoic orogens, such as the Tibetan Plateau, generally have weak lithospheric strength, while the Precambrian orogens have high strength (Deng and Tesauro, 2016;Sun et al., 2013). The results are also supported by the S-wave velocity structure across the margins of the Tibetan Plateau and adjacent blocks (Tian et al., 2021). This is not only related to the cooling of the lithosphere but also to the composition of the continental lithosphere and the distribution of radioactive elements (Mareschal and Jaupart, 2013;Sun et al., 2022). In this study, for the reference model (case R), the initial lithospheric thickness of the craton does not change. The parameters of case R are given in Table 1. When we consider the effect of the lithospheric strength heterogeneities on the strain partitioning, a weak orogen is inserted in the craton. According to the nature of the orogens in the periphery of the craton or the interior of the craton, the lithospheric strength or thickness contrasts are considered. This difference is inherited from older periods, such as the huge difference in lithospheric thickness of the North China craton after the Mesozoic (Zhang et al., 2019). In addition, the effect of crustal strength contrast is considered. Rheological parameters adopted in these numerical models follow in Table 1. 3 Results

Reference model (case R)
For case R, the Tibetan Plateau and craton correspond to high and low lithospheric viscosity, respectively ( Figure 2B). With continuous convergence (20 mm/yr on the left side of the model), the results show that the lithosphere of the Tibetan Plateau has undergone significant deformation, and a high strain rate is localized on it ( Figures 3A-D). The mantle lithosphere of the Tibetan Plateau is delaminated or subducted after encountering the obstruction of the craton (Figures 3C,D). The craton does not deform significantly, which is indicated by a low strain rate or low convergence velocity ( Figure 3E). The crustal thickening and high strain rate are mainly distributed on the side of the Tibetan Plateau near the boundary between the plateau and the craton. The free-slip boundary on the right side of the model means that the length of the craton could extend to infinity. Under this high rheological contrast ( Figure 2B; Figure 3D), these results indicate that the expansion of the Tibetan Plateau is blocked by the craton. If the craton is uniform and infinite, the expanding plateau will not have a far-field effect under the India-Eurasian collision.

Effect of rheological heterogeneity of the lithosphere (group A)
To further study the effect of lateral rheological heterogeneity on lithospheric deformation and strain partitioning. Three more models (cases A, A1, A2) were established ( Figure 4) for comparison with case R. In these cases, a 150 km wide orogen is inserted in the craton. In case A, the weak orogen was imposed the same rheological parameters as the Tibetan Plateau (Table 1). The rheological parameters of case A1 are similar to case A, except that the orogen was given a lower plastic yield strength (10 MPa). For case A2, the mantle lithosphere of the Tibetan Plateau is given a higher strength (Dry olivine replaces wet olivine according to Table 1)

Frontiers in Earth Science frontiersin.org 05
Furthermore, if the lithospheric strength of the orogen is decreased, the deformation is more likely to be transmitted to the orogen outside the craton (case A1 in Figures 4C,D). In case A2, the strong mantle lithosphere of the Tibetan Plateau collides with that of the craton, which causes a high strain rate and delamination of the mantle lithosphere under the orogen (Figures 4E,F).
According to the results of the viscous thin-shell model (England and Houseman, 1986), during the outward expansion of the Tibetan Plateau, strong deformation and high strain rates are localized on the weak Tibetan Plateau. For a homogeneous crust, high strain gradually decreases away from the collision boundary and is hard to transfer to the orogen outside the craton. However, the lateral strength heterogeneity of the lithosphere changes this feature, making it easier for the high strain to concentrate on the Tibetan Plateau or orogens with weak strength (England and Houseman, 1985). In our two-dimensional profile model, such conclusions also apply to the horizontal deformation characteristics and strain distribution of the lithosphere. In addition, the deformation of the lithosphere also affects the motion of the asthenosphere. For example, the deformation of the mantle lithosphere will cause small-scale mantle convection below the orogen ( Figure 4F).
We compared the surface topography and strain rate of these three models (group A) with the observations across the northern margin of the Tibetan Plateau (AA' profile in Figure 1B). The comparisons ( Figure 5) show that the results of group A are generally consistent with the actual deformation feature. The west Kunlun Shan and Tian Shan with high landforms are located on the south and north sides of the Tarim basin, which is similar in the high topography features distributed over both sides of the craton obtained by the numerical models ( Figure 5A). Especially for case A2, the steep topography feature is fits well with the observation profile. High values of the maximum shear strain rate ( Figure 5B) and dilatation strain rate ( Figure 5C) from the observations (Kreemer et al., 2014;Wang and Shen, 2020), show that Tian Shan are under strong compressive deformation. These are generally consistent with the results of numerical models. The features of observed topography and strain rate of the Than Shan are in better agreement with the results of case A1 or case A2. This implies that the lithospheric thickening may be taking place beneath the Tian Shan. However, the dilatation strain rate calculation shows that the Tibetan Plateau is under a compression state, which is not consistent with the observed extension state ( Figure 5C). This may be related to the shallow crustal extension caused by gravity collapse for high topography (Liu Frontiers in Earth Science frontiersin.org 06 and Yang, 2003), because the current observed strain rate is calculated from the GPS observation on the surface. It represents the threedimensional deformation characteristics of the shallow crust.

Effect of lithospheric thickness on deformation (group B)
Given the same rheological parameters, the initial thickness of the lithosphere also affects the strength of the lithosphere because of the difference in thermal structure. For instance, the craton and its internal orogen have different strengths because of the varied initial thicknesses of the lithosphere. This is related to the evolution of craton or orogens. Therefore, it is necessary to discuss the influence of the initial structure of the lithosphere on the deformation and strain distribution.
Based on case R, we set up two comparable cases (case B1 and case B2) by changing the initial lithospheric thickness of the craton ( Figure 6). For case B1, the initial lithospheric thickness of the craton is not constant. The initial lithospheric thickness in the left part of the craton is 130 km, which is thicker than that of the right part (100 km) (Figures 6C,D). For case B2, a much thinner initial lithospheric thickness (80 km) of the right part of the craton is imposed (Figures 6E,F). The results show that the difference in the initial thickness of the lithosphere, which affects the rheological strength (Supplementary Figure S2), also has an important influence on the deformation and strain partitioning. In case R, the craton with the same initial lithospheric thickness has higher lithospheric strength than the Tibetan Plateau. During the expansion of the Tibetan Plateau, the deformation and high strain are mainly concentrated in the Tibetan Plateau. With the initial lithospheric thickness of the right part of the craton decreasing, the effective viscosity also decreases (Supplementary Figure S2). The strain rate localized on the right part of the craton increases accordingly  Comparisons of the results of numerical models (group A) with the observed surface topography (A), the maximum shear strain rate (B), and the dilatation strain rate (C) across the northern margin of the Tibetan Plateau (AA' profile in Figure 1B). Gray areas in (A) indicate the observed topographic change. Gray dots in (B-C) indicate the observed strain rates (Kreemer et al., 2014;Wang and Shen, 2020).

Frontiers in Earth Science
frontiersin.org ( Figures 6C-F). In case B2, if the initial lithospheric thickness of the right part of the craton is reduced by 50 km compared to case R, the strain rate localized on the right part of the craton is similar to that of the Tibetan Plateau ( Figures 6E,F). The strain rate localized on the left part of the craton is similar to that of case R. The surface topography and strain rate of these three cases (cases R, B1, and B2) are compared with the observations across the eastern margins of the Tibetan Plateau (BB' profile in Figure 1C). The results of the reference model (case R) are consistent with the observations of this profile. Crossing the boundary between the plateau and the craton, the topography and strain rate changed abruptly (Figure 7). When the lithosphere thickness of the right part of the craton decreases (case B1, B2), the topography and strain rate of the right part of the craton will increase ( Figures 7A,B). High topography and the maximum shear strain rate correspond to the low effective viscosity of the lithosphere (Supplementary Figure S2), which is determined by the initial lithosphere thickness. Our model shows that when the initial lithosphere thickness of the right part of the craton is 100 km (case B1), the calculated surface topography and the maximum strain rate agree with the observation results ( Figures 7A,B). This is consistent with the thermal lithospheric thickness ( Figure 1C). However, the comparison of dilatation strain rate shows the compressive state from the simulation is stronger than that of observation ( Figure 7C).

Effect of crustal rheological contrast (group C)
The above cases include the overall rheological contrast of the lithosphere, but the rheological contrast of the crust has not been considered in detail. In this group of cases, the effect of the crustal rheological contrast on the plateau expansion and far-field effect is considered. In case C ( Figures 8A,B), we imposed a weak crust for the craton (Supplementary Figure S3), whose parameters are the  Comparisons of the numerical models (cases R, B1, B2) with the observed surface topography (A), the maximum shear strain rate (B), and the dilatation strain rate (C) across the eastern margins of the Tibetan Plateau (BB' profile in Figure 1C), respectively. Gray areas in (A) indicate the observed topographic change. Gray dots in (B-C) indicate the observed strain rates. The data of topography are consistent with those shown in Figure 1. Strain rate data is from (Kreemer et al., 2014;Wang and Shen, 2020 Figure S3B), the high strain rate can be extended to the interior of the craton under the expansion of the Tibetan Plateau ( Figure 8B). The strain rate results show that the deformation process is gradually extended from the plateau and craton boundary to the craton interior. As shown in group A, if there is a weak orogen in the craton ( Figures 8C,D), the deformation of the orogenic belt is relatively intense and a high strain rate is localized on it (case C1). Furthermore, if a thin initial lithosphere is imposed for the right part of the craton (case C2), it is conducive to the high strain rate localized on it. In the upper mantle, edge-driven convection will be caused by the large lithospheric thickness of the craton (Figures 8E,F). The surface topography and strain rate of these three cases (cases C, C1, and C2) are compared with the observations across the northeastern margins of the Tibetan Plateau (CC' profile in Figure 1D). For the northeastern margin of the Tibetan Plateau, the smooth surface topography and the diffuse strain rate can be generally fitted only when a weak crustal strength of the craton is given (Figure 9). If the orogen within the craton is given low lithospheric strength, then the surface topography and strain rate are consistent with the characteristics of the Trans-North China orogen (case C1 in Figure 9). This is reasonable because the composition between the ancient orogenic belt and stable craton are still different. However, different compositions may also produce the same rheological strength (group B). In our current model, the factor affecting the results can be attributed to the rheological contrast (Supplementary Figure  S3). In case C2, under the convergent of the Tibetan Plateau, low viscosity and large difference of lithospheric thickness leads to strong deformation of the orogen in the craton. High topography and sharp change of the strain rate crossing the orogen in the craton is not consistent with the observed strain rate (case C2 in Figure 9). Thus, we infer that the lithospheric

FIGURE 9
Comparisons of the numerical models (cases C, C1, C2) with the observed surface topography (A), the maximum shear strain rate (B), and the dilatation strain rate (C) across the northeastern margins of the Tibetan Plateau (CC' profile in Figure 1D), respectively. Gray areas in (A) indicate the observed topographic change. Gray dots in (B-C) indicate the observed strain rates. The data of topography and strain rate are consistent with those shown in Figure 5.

Frontiers in Earth Science
frontiersin.org deformation of the Trans-North China is not strong as that shown in case C2.

Implication for the far-field effect of the India-Eurasian collision
The Tibetan Plateau is blocked by cratons during its expansion, but whether it has far-field effects on orogens outside or within the cratons has been debated. Using numerical models in this study, we considered the effects of strength heterogeneity of the lithosphere and contrast in lithospheric thickness on the deformation characteristics and strain distribution during the expansion of the plateau. These results can help us understand the far-field effects of the India-Eurasian collision.
Across the northern margin of the Tibetan Plateau and the Tian Shan orogen (profile AA' in Figure 1B), the lithospheric deformation can be understood by the lateral heterogeneity of lithospheric strength. Our numerical results show that lateral heterogeneity of lithospheric strength affects greatly lithospheric deformation patterns and strain partitioning (Figures 4, 5). When the expanding Tibetan Plateau is blocked by the Tarim, the hard Tarim can transmit the strain to Tian Shan orogen (Figure 4). These results are consistent with previous thin-shell model studies (England and Houseman, 1985;Neil and Houseman, 1997). The pre-existing weak Tian Shan orogen was remotely influenced by the India-Eurasian collision (Kong et al., 1997). Especially when the subducted Indian mantle lithosphere collided with the mantle lithosphere of the Tarim (Li et al., 2008;, the uplift of the Tian Shan orogen was more intense. For case A2, if the weak mantle of the Tibetan Plateau was replaced by the subducted Indian mantle lithosphere ( Figures 4E,F), lithospheric thickening occurred strongly on the Tian Shan orogen. According to the observation of surface strain rate and geomorphological characteristics ( Figure 5), our numerical results show that the lithospheric thickening beneath Tian Shan orogen is consistent with the characteristics observed in some geophysical observations (Chen et al., 1997;He and Santosh, 2018). Therefore, the late Cenozoic reactivation of the Tian Shan orogen was controlled by the India-Eurasian collision, especially the collision of the mantle lithosphere between India and Tarim since the late Cenozoic (Hendrix et al., 1994;Qin et al., 2022). But the low lithospheric strength of the Tian Shan orogen is the premise that it can be remotely affected by the India-Eurasian collision.
However, when the lithospheric strength of the orogenic belt in the craton is relatively high, the India-Eurasian collision is difficult to affect (Figures 6, 7). This result can be used to understand the convergence deformation between the Tibetan Plateau and South China block (profile BB' in Figure 1C). Within the South China block, there is the Jiangnan orogen formed by the splicing of the Yangtze craton and the Cathaysian block in the Neoproterozoic (Li, 1999;Li et al., 1995). Since these orogens was formed during the collision of two blocks, its strength is similar to that of the Yangtze craton and the Cathaysian block. Furthermore, its strength is higher than the orogen formed in the Paleozoic, such as Tian Shan orogen (Allen et al., 1993). When such a high strength ancient orogenic belt exists in the craton, the strain caused by India-Eurasian collision cannot be transited to it (Figures 6B,C; Figures 7A,B). That means it is difficult for this ancient orogenic belt to be reactivated like the Tian Shan orogen. In particular, if the craton and its interior orogenic belt remain thick lithosphere, the Tibetan Plateau will be severely constrained, making it difficult to expand further ( Figures  6A-D). Deformation and high strain rate are mainly localized on the Tibetan Plateau, and the mantle lithosphere of the extended Tibetan Plateau will undergo delamination, which is very similar to the deformation characteristics of the eastern margin of the Tibetan Plateau (Liu et al., 2014;Wang et al., 2014b). Our numerical model results suggest that the lithospheric thickness also affects strain distribution, although not as strong as the preexisting weak lithospheric strength Figures 7B,C). The current thickness of the mantle lithosphere in the Cathaysia block in the east of the South China block has decreased significantly (An and Shi, 2006;Gao et al., 2022;Li et al., 2013). Although thin lithospheric thickness contributes to high strain rate distribution, our results show that the high lithospheric strength of the craton will still restrict the expansion of the Tibetan Plateau and its far-field effect.
For the northeastern margin of the Tibetan Plateau (profile CC' in Figure 1D), more complex lithospheric rheology is needed to fit the observed topography and strain rate. Crossing this profile, the smooth surface topography and the diffuse strain rate require low crustal rheological contrast between the Tibetan Plateau and the North China craton (Figures 8, 9). In addition, our results show that the formation of the smooth surface topography needs a weak mantle lithosphere of the west part of the North China craton (or the Ordos block). Although the Trans-North China orogen was formed by the splicing of the eastern and western blocks of the North China craton in the Paleoproterozoic (Zhao, 2001;Zhao et al., 2000;Zhao et al., 2005), its strength should be as high as that of the Jiangnan orogen. However, our numerical results show that the high topography of the Trans-North China orogen needs a weak lithospheric strength (Figures 9B,C). This may imply that the Trans-North China orogen had undergone significantly transformation during the Mesozoic (Zhu and Xu, 2019). For example, during the Mesozoic, the mantle lithosphere in the eastern North China craton was significantly thinned from about 120 km to the current 80 km (Chen et al., 2009;Zhang et al., 2014), and the strength of the Trans-North China orogen was reduced. In addition, differences in lithosphere thickness can cause the edge-driven convention beneath the North China craton (King and Anderson, 1998;King and Ritsema, 2000;Sun and Liu, 2023). The superposition of edge-driven convention and asthenosphere flow caused by the expansion of the Tibetan Plateau will further reduce the strength of the lithosphere beneath the Trans-North China orogen ( Figures  8B,D,F). Thus, the lateral growth of the Tibetan Plateau is most likely to occur in the northeast margin and may affect the Trans-North China orogen.

Frontiers in Earth Science
frontiersin.org

Uncertainties of the model
In this study, the push boundary condition on the left side of the model presents the convergence of the Indian and Eurasia plates. However, the convergence of the Eurasia-Pacific plates is also an important driving force that cannot be ignored (Northrup et al., 1995). The subduction of the Paleo-Pacific (Izanagi) plate is regarded as the dominant factor that induced lithospheric thinning of the North China craton in the Mesozoic (Liu et al., 2021a;Zhu and Xu, 2019). The subsequent subduction and roll-back of the Pacific plate lasted to the present and trigged back-arc extension in east China (Liu et al., 2017a). Especially, the Cenozoic big mantle wedge caused by the stagnation of the subducting Pacific slab has contributed to the intraplate volcanism and rifting in East Asia (Liu et al., 2017b;. The spatial patterns of the fast polarization directions of vertically traveling shear waves also indicate complex upper mantle deformation in east China (Zhao and Xue, 2010). Strain rates of group C models in this study show the deformation of the right part of the craton is under compression, which is not consistent with the observed extensive strain rates in the east part of the North China craton ( Figures 8G,H). This may be related to the fact that the pulling force caused by the roll-back of the Pacific plate is not considered in our model (Liu et al., 2021b;Yang et al., 2018). If a tensile force is provided at the right boundary in the numerical model, the lithosphere of the surrounding block will be under extension state according to recent study (Sun and Liu, 2023).
At the same time, there are large differences in the lithosphere thickness between the east and west of the North China craton after being reconstructed. Deep geophysical results indicate that crossing the Ordos, there are low-velocity anomalies and low electrical resistivity under the Trans-North China orogen (Guo et al., 2016;Lei, 2012;Li et al., 2018;Yin et al., 2017), indicating upwelling and downwelling of the asthenosphere beneath the North China craton. Although the existence of the Philippine plate blocks the impact of the westward subduction of the Pacific plate on the South China block (Liu et al., 2021b), similar tomography results are also observed in the upper mantle beneath the South China block (Gao et al., 2022;Ward et al., 2021;Wei et al., 2012). This may imply the existence of smallscale convection in the upper mantle of east China because of the difference in lithospheric structure Anderson, 1995, 1998). These studies show that the complicated lithospheric and upper mantle deformation of east China may be the interaction of the subduction of the Pacific plate, the India-Eurasian collision and small-scale convection (Sun and Liu, 2023).
Furthermore, the two-dimensional model here cannot consider the lateral deformation or mantle flow caused by the India-Eurasian collision. On the northeastern margin of the Tibetan Plateau, the extended plateau not only converges with the western margin of the North China craton but also experienced significant strike-slip and developed a series of pull-apart basins (Zheng et al., 2013). In the east of the Tibetan Plateau, current GPS and geophysical observations show that the plateau material is extruded along the southeastern margin of the plateau and the Qingling orogen (Guo and Chen, 2017;Wang and Shen, 2020). The velocity boundary conditions (20 mm/yr) we impose on all models are the same. Actually, for the region of the eastern and northern margins of the Tibetan Plateau, only the velocity perpendicular to the block boundary should be considered in the two-dimensional model. In this way, the convergence velocity is lower than the current given velocity boundary conditions. Thus, the strain localized on the orogens in the craton may be smaller than present. These lateral deformations may reduce the far-field effects of the India-Eurasian collision because it does not need to overcome great resistance in the case of the pre-existing weak zone (Chen et al., 2020;Kong et al., 1997). Therefore, the Trans-North China orogen and Jiangnan orogen would be more stable in that case. The results of this study will not be changed.
Finally, to test the model in the future, systematic comparisons between the model results and observations should be strengthened. Besides topography and strain rate, other surface geology and geophysical data from multiple disciplines, such as the crustal-lithospheric structure and thermal structure, should be used as much as possible. However, there are still two problems in the comparison of lithosphere structure. On the one hand, the lithosphere structure obtained by varied methods is quite different, such as the lithosphere thickness obtained by temperature and seismology can vary greatly [An and Shi., 2006;Pasyanos et al., 2014]. On the other hand, the lithosphere structure obtained by simulation is highly dependent on the initial structure, which will cause great uncertainty in comparison.

Conclusion
By constructing two-dimensional thermal-mechanical numerical models, this study investigated the effects of the rheological heterogeneities and thicknesses of the lithosphere on the deformation of the Tibetan Plateau and neighbouring regions. These results help us understand the far-field effects of the India-Eurasian collision. The main conclusions are as follows.
(1) The lateral rheological heterogeneities of the lithosphere have an important influence on the deformation and strain partitioning during the lateral expansion of the Tibetan Plateau. When the expanding Tibetan Plateau is blocked by a strong peripheral craton, if there is a weak orogen outside the craton, the craton will transmit the strain to the peripheral orogen and reactivate it, such as the Tian Shan orogen. The collision of the subducted Indian mantle lithosphere with the Tarim mantle lithosphere would make the uplift of the Tian Shan orogen more intense.
(2) If there is an ancient orogenic belt in the peripheral craton, the deformation is mainly concentrated in the weak-strength Tibetan Plateau, and the lateral expansion of the plateau is limited by the craton. For example, the Jiangnan Orogen within the South China block is hard to be reactivated by the collision between the Indian and Eurasia plates.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.