Lower Crustal Rheology Controls the Development of Large Offset Strike‐Slip Faults During the Himalayan‐Tibetan Orogeny

The mechanism of crustal deformation and the development of large offset strike‐slip faults during continental collision, such as the India‐Eurasia zone, remains poorly understood. Previous mechanical models were simplified which are either (quasi‐)2‐D approximations or made the a priori assumption that the rheology of the lithosphere was either purely viscous (distributed deformation) or purely localized. Here we present three‐dimensional visco‐elasto‐plastic thermo‐mechanical simulations, which can produce both distributed and highly localized deformation, to investigate crustal deformation during continental indentation. Our results show that large‐scale shear zones develop as a result of frictional plasticity, which have many similarities with observed shear zones. Yet localized deformation requires both a strong upper crust (>1022 Pa·s) and a moderately weak middle/lower crust (~1020 Pa·s) in Tibet. The brittle shear zones in our models develop low viscosity zones directly beneath them, consistent with geological observations of exhumed faults, and geophysical observations across active faults.


Introduction
The Himalayan-Tibetan orogen is a spectacular example of continent-continent collision, which resulted in pervasive crustal deformation in Asia (Yin & Harrison, 2000). Such deformation includes several large-scale, extending 1,000 km long, strike-slip faults (Figure 1) that developed through the collision history (Molnar & Tapponnier, 1975;Tapponnier et al., 1986). The mechanisms of fault development and the partitioning of crustal deformation, however, are poorly understood. Several endmember models have been proposed for lithospheric deformation in Tibet and surrounding regions, including an extrusion model suggesting that deformation is primarily localized along large-scale strike-slip faults and thrust zones (Replumaz & Tapponnier, 2003;Tapponnier et al., 1982Tapponnier et al., , 2001; thin viscous sheet models suggesting that shortening is distributed and vertically coherent within the lithosphere (England & McKenzie, 1982;Houseman & England, 1986); and channel flow models implying that a thickened weak channel develops in the lower crust which spreads the deformation front outward from the plateau (Clark & Royden, 2000;Royden et al., 1997).  Gan et al. (2007); orange shaded regions are low seismic velocity zones (Bao et al., 2015;Liu et al., 2014;Wittlinger et al., 1998); purple shaded regions are high electrical conductivity zones (Bai et al., 2010;Dong et al., 2016;Le Pape et al., 2012;Wei et al., 2014;Zhang et al., 2015). The black dashed zone is roughly our model region.

Geophysical Research Letters
Deformation style varies from normal faulting to thrust and strike-slip faulting across the plateau and its periphery, indicating complex deformation processes during the India-Eurasia collision. Present-day strain partitioning is locally mainly localized along major strike-slip faults, supporting an extrusion model . Geophysical studies show that at several locations low seismic velocity zones are coincident at depth with high electrical conductivity zones around the large offset strike-slip faults (Agius & Lebedev, 2014;Bai et al., 2010;Baumann & Kaus, 2015;Nelson et al., 1996;Wu et al., 2019) (Figure 1). These anomalous zones in the middle-lower crust are generally attributed to high temperature, the presence of aqueous fluids, and/or partial melting Liu & Hasterok, 2016), which will all significantly weaken the crustal strength (Kaus, 2016;Liu & Hasterok, 2016) beneath the faults. Furthermore, the present-day deformation mode in the lower crust remains a matter of debate: Copley et al. (2011) argue in favor of a strong lower crust beneath Tibet to explain strike-slip faulting in the north and normal faulting in the south, which contrasts with the channel flow model (e.g., Clark & Royden, 2000). As geophysical data give a snapshot of present-day crustal and mantle structures, the properties before collision remain unclear, and its role in the evolution of plateau growth and crustal deformation remains ambiguous.
The development of crustal deformation and topography are strongly dependent on the rheology of the crust, which is the fundamental distinction among the endmember models listed above. Nevertheless, most previous models were either fully plastic for the block model or fully viscous for thin sheet and lower crustal flow models. Recent numerical models mostly focus on the development of topography (Chen et al., 2017;Lechmann et al., 2011;Pusok & Kaus, 2015) and subduction dynamics (Capitanio, 2020;Li et al., 2013;Replumaz et al., 2016;Sternai et al., 2014), emphasizing the effect of variable lateral and vertical lithosphere strength across multiple terranes (Bischoff & Flesch, 2018;Chen et al., 2017;Huangfu et al., 2018). Such models provide valuable insights in the dynamics of the lithospheric deformation, but the limited resolution of these models inhibits resolving strongly localized deformation within the crust in a three-dimensional (3-D) framework.
Here, we employ high-resolution three-dimensional thermo-mechanical numerical models to investigate the role of crustal rheology on the formation of large strike-slip fault during the indentation of India into Asia.

Model Description
We employ the 3-D parallel thermo-mechanical numerical code LaMEM , which is based on a finite difference staggered grid discretization combined with a particle-in-cell technique using a conservative velocity interpolation approach (Pusok et al., 2017). It solves the energy, momentum, and mass conservation using a geometric multigrid approach (see supporting information for detail). Our simulations are performed for a 1,500 km × 2,400 km × 55 km ( Figure 1) region (using 512 × 768 × 32 cells) and employ a simplified tectonic structure with nonlinear visco-elasto-plastic rheologies (Pusok et al., 2018;Wang et al., 2019), designed to investigate crustal deformation of an India-like indenter into Asia. The 40-kmthick crust is composed of a 20-km upper and a 20-km-thick lower crust in Tibetan Block (TB) region, respectively. The vertical grid resolution has been increased to 1.25 km in the upper crust to capture localized deformation (Duretz et al., 2020). The mechanical boundary conditions are free slip for all boundaries, but to mimic a free surface condition in finite difference scheme, we use a 15-km "sticky-air" layer as well as with a free surface stabilization algorithm (Crameri et al., 2012;Kaus et al., 2010). The strong indenter collides with a constant convergence (5 cm/yr) into TB ( Figure 1). Three weak seeds are placed in the TB region (the effect of their locations is tested, as described below). For simplicity, we employ an initial linear temperature profile with a fixed temperature of 0°C at the surface and 600°C at the bottom, while other boundaries have zero thermal flux.
The rheology of rocks (Table S1) is assumed to be visco-elasto-plastic, and the total deviatoric strain rate is given by where _ ε vis ij , _ ε el ij , _ ε pl ij denote viscous, elastic, and plastic strain rate, respectively. η eff is the effective viscosity, G elastic shear modulus, t time, _ λ plastic multiplier, Q plastic flow potential, and σ ij = − P+τ ij the total stress. The power law viscosity is given by 10.1029/2020GL089435

Geophysical Research Letters
where f is the scaling factor used in the simulations (Lu et al., 2016), B N is the prefactor, n is the stress is the second invariant of the viscous strain rate tensor, E, V are activation energy and activation volume, respectively, and R is the gas constant. We assume Drucker-Prager yield criteria with which rocks fail plastically if the second invariant of the deviatoric stress (τ II ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 0:5τ ij τ ij p ) exceed the yield stress τ y (i.e., τ II > τ y ). For Drucker-Prager non-associated plasticity this is expressed as where P is the pressure, C is the cohesion, and ϕ is the effective friction angle. Plastic strain weakening of the friction angle is applied in the model which linearly decreases the friction angle from 30°to 3°b etween a strain of 0.5 and 1.5 (Table S1).

Results
We use a simplified tectonic geometry which includes the TB and three strong blocks in the north (SBN), the south (SBS, the Indian indenter), and the east (SBE), respectively ( Figure 1). We focus on the formation of large-scale faults during continental collision, and the influence of crustal strength is discussed first, after

Geophysical Research Letters
which we present additional experiments to test the influence of shear heating and the location of the weak seeds. We start with explaining the reference model (M3) which best matches the observations.

Crustal Rheology Controls Strike-Slip Fault Formation
Upon the indentation of strong block SBS toward north (y direction), an intense simple shearing between the indenter and the adjacent TB block induces the first shear zone (or fault)-SZ1 (Figure 2). Six more shear zones successively develop within the TB block at later stages along either preset weak seeds or the boundary of strong blocks. The deformation is strongly concentrated along those narrow shear zones in the upper

Geophysical Research Letters
crust, while it is slightly more diffuse in the lower crust. The shear zones show a complex deformation pattern, exhibiting predominantly strike-slip motion, but also small amounts of thrust and/or normal slip within (e.g., SZ4, SZ5, and SZ7 in Figure 2). In front of the indenter, the lower crust is thickened and high topography develops, while the upper crust is extruded in east-and southeast-ward directions, rotating around a stationary point ( Figure S1). Both thrust and normal faults occur in the frontal and side edges, respectively, due to compression and gravitational collapse. Escaping crust is blocked by SZ5, leading to thrust faulting in the northeast. When the indenter moves further northward, SZ5 rotates away from the indenter, extruded crustal material which is subsequently bounded by SZ7 formed along SBN. The shear zones are long, some being 1,000 km long, either continuous or segmented (Figure 3). They are either mostly left-lateral (SZ3, SZ5, and SZ7) or right-lateral (SZ1, SZ4, and SZ6). Slip displacements are large: SZ1 has offsets of~100-400 km, SZ3~133 km, SZ5~382 km, and SZ7~378 km. There are also small widespread faults in the southeastern edge, primarily limited to the shallow upper crust at the early collision stage. In the southeast, no large shear zones develop within the crustal block between SZ1 and SZ3.
In order to test the impact of crustal strength on the dynamics, we perform four endmember simulations in which we test the effect of a very weak lower (M5; f uc = 1, f lc = 0.01, f being the scaling factor of viscosity in Equation 2) or a weak upper crust (M12; f uc = 0.01, f lc = 1), a strong lower crust (M1; f uc = 1, f lc = 100), and a case with a weak middle crust and a moderate lower crust (3Ly; f mc = 1) (see Figure 3).
The simulation with a very weak lower crust (M5) shows a distinct behavior from others in that the upper crust is strongly thickened in front of the indenter, and the lower crust is almost exhumed to the surface beneath SZ5. SZ3 is not developed in Model M5, while SZ2, which is absent in the other models, is well developed. On the other hand, the case with a strong lower crust (M1) has far less remarkable shear zones, and a more pronounced thickening of the lower crust, which expands toward the east and southeast. Yet there is no correlation between the development of shear zones and thickening of the lower crust in Model M1. This also occurs for the case with a very weak upper crust (M12) in which shear zones are poorly developed, even though the lower crust undulation is similar as in Model M3. The case with a weak middle crust (3Ly) showed roughly similar shear zones as the reference simulation, but with less lower crustal elevation beneath the shear zones.

Initiation of Shear Zones
In the reference model, shear zones develop either along strength contrasts of two terranes or from local heterogeneities such as the initially imposed weak seeds. To have a better insight into the development of the shear zone, we describe SZ5 as an example (Figure 4). As a result of strike-slip motion, the lower crust moves upward and leads to~6 km uplift of the lower crust, contemporaneously with the localization of SZ5. SZ5 initiates at~12-km depth around 19 Myr and further localized at~19 km. The whole crust is strongly deformed, resulting in a drastic undulation of the upper-lower crust interface, particularly in regions where shear localization occurs. The localization initiation depth can be either in the uppermost crust (<10 km) or in the middle crust (~15 km), but only those faults that are in the vicinity of large viscous terranes or that initiate from weak seeds extend to the lower crust depth and develop large-scale shear zones. This thus suggests that stress or mechanical heterogeneities play a key role in the initiation of such fault zones.
The model without shear heating (M13) exhibits less pronounced shear zones which generally locates around the strong terranes (Figure 5a). This points to a critical influence of shear heating (Kiss  Leloup et al., 1999;Popov et al., 2012) and strong terranes on the strain localization. If we remove the strong block SBE (M14 in Figure 5), SZ3 and SZ4 are not generated and SZ5 is very poorly developed. To test the influence of the weak seeds, four additional experiments are performed: We either remove the weak seed(s) or place it at a different position. The model results show that strike-slip faults still form without preset weak seeds, but they are shorter, as it is the case for SZ5 in the reference model, indicating a relatively weak correlation of the weak seed 1 with SZ5.

Discussion and Conclusions
Although the model geometry is simplified, these models are useful to understand the dynamics of crustal deformation during Himalayan-Tibetan orogeny.  Figure S1). A short indenter allows the crust to flow to the west after its passage, which results in an early developed SZ1 while the Sagaing Fault which might be affected by retreating Indian Plate was initiated at~22-16 Ma (Searle et al., 2007, and reference therein). Interestingly, SZ2 that does not appear in the reference model, yet is developed when the lower crust is even weaker (<10 19 Pa·s), indicating that the formation of Red River Fault is probably related to a locally weak crust. In contrast, deformation is less localized in the upper crust with no brittle faults if the upper crust is too weak. This implies that a seismogenic layer should be sufficiently strong. However, if the lower crust is too strong (>10 22 Pa·s), the whole crust undergoes pure shear thickening which results in coherent deformation and less developed shear zones. Thus, the absent SZ2 may be potentially explained by different deformation styles due to along-strike rheological heterogeneities between the east-west (Bischoff & Flesch, 2019;Chen et al., 2017;Zheng et al., 2020) and the north-south (Huangfu et al., 2018) Tibetan terranes. Notably, the effective viscosities of upper and lower crust in our reference model are consistent with previous models which account for normal faulting and viscous buckling (Bischoff & Flesch, 2018;Flesch et al., 2001) in Tibet.
A comparison between Models M3 and M1 (3Ly) shows that a weak lower (middle) crust facilitates shear zone development (Figure 3). The locations of shear zones have a positive relationship with lower crust upwelling. The shear zones are usually localized at 15-to 20-km depth, with shear heating being the major heat source that lowers in situ strength and further promotes strain localization (Kaus & Podladchikov, 2006;Kiss et al., 2019;Leloup et al., 1999), while the initial thermal structure of the crust has a limited influence on the localization ( Figures S2 and S3). The lower crust is uplifted beneath the shear zones, which is a robust feature in our simulations, further increasing the local temperature and facilitating strain localization. A similar increase in temperature and uplift of the middle crust has been documented in exhumed roots of large strike-slip faults (Leloup et al., 1999, and references therein). The uplift of the crustal flow probably indicates that the long-lived shear zones are predominantly transtensional (Le Pourhiet et al., 2014) or transpressional (Popov et al., 2012) during their formation. Our model results show that the shear zone close to a strong block likely exhibits reverse fault, transpression, and strike-slip fault; otherwise, it shows transtension and strike slip ( Figure S4). Therefore, our model results predict a thickened lower crust with low seismic velocity beneath the major shear zones. This is consistent with seismic tomography and magnetotelluric observations of such low seismic velocity (Agius & Lebedev, 2014;Bao et al., 2015;Jiang et al., 2014;Liu et al., 2014;Wang et al., 2018) and high electrical conductivity zones (Bai et al., 2010;Wei et al., 2001) in the middle/lower crust beneath the Xianshuihe and Kunlun faults. Similar features have also locally been observed along the Altyn-Tagh Fault (Wittlinger et al., 1998;Zhang et al., 2015). Two remarkable seismic low-velocity channels (Figure 1) in southeast Tibet correlate with crustal thickening. It is worth noting that most large earthquakes occur within the boundaries of these two channels (Bao et al., 2015), which corroborates our finding that the formation of large-scale shear zones require/promote a weak and flowing middle/lower crust. Indeed, other numerical models also show that a weak middle to lower crust is required for the development of the N-S rifting systems in southern Tibet (Pang et al., 2018), as well as for the development of a wide plateau (Chen et al., 2017).
The reactivation of inherited weakness, such as ancient suture zones, is often proposed to control the location of lithosphere deformation during continental collision Sun et al., 2018). However, few of the large-scale strike-slip faults in Asia geographically follow previous suture zones, with the notable exception of the Kunlun Fault (Tapponnier et al., 1986) (Figure 1). Previous models use weak zones that cannot sustain even small stresses (probably << few tens MPa) and thus localize deformation along the whole weak belt simultaneously, which is inconsistent with the observed diachronous initiation along strike-slip faults, such as Kunlun Fault which initiates in the west first and then propagates to the east (Duvall et al., 2013). In contrast, the presence of craton-like terranes in the periphery of Tibet, which sustains large stresses upon collision (Lu et al., 2011), does play a dominant role for the initiation of large-scale faults in our model (Figure 5b), consistent with geological observations in Tibet region.
The surface dynamics of continental extrusion during plate convergence is usually attributed to the mantle flow (Capitanio, 2014;Capitanio & Replumaz, 2013;Sternai et al., 2014). This would postulate that the mantle traction force is sufficiently large (Capitanio, 2020) and/or that the crust and mantle are coupled (Sol et al., 2007;Wang et al., 2008). However, our models show that a strong lower crust prohibits localized deformation, suggesting that the middle/lower crust is at least locally weak beneath the faults.
In conclusion, our simulations show that a moderately weak middle/lower crust is required to develop long-lived large-scale strike-slip faults as observed in Asia. The ongoing continental indentation pushes crustal blocks to "escape" toward northeast and southeast Tibet, causing tectonic forces and localized deformation along strike-slip faults that accommodate displacements. Deformation is primarily localized in shear zones which cut through the entire crust, whereas less pronounced deformation is localized in small faults in the rigid upper crust.

Data Availability Statement
The numerical source code used in this study is an open source online and also available from the link https://figshare.com/articles/LaMEM_TibetCollision/12571955, and the model setup to reproduce the results can be found on https://figshare.com/articles/inputdataTibet/12573236 or request from the corresponding author.