Analytical and numerical simulations of uplift processes at the Tibet-Sichuan boundary

Previous studies have shown that the uplift of Tibetan plateau started in response to the collision of Indian plate and Eurasian plate. During this process, the crust of Tibetan plateau has been greatly thickened which leads to significant elevations. The elevation gradient is extremely large at the east boundary of Tibetan plateau where Longmenshan fault exists, dropping from 4500 to 500 m within a distance of 100 km, while it is more gentle at the south and north sides of Sichuan basin. Such a difference of elevation gradient has been explained with a crustal channel flow model. However, previous crustal flow models consider the thickness of the lower crust as a constant which is highly simplified. Therefore, it is essential to build a more realistic crustal flow model, in which the thickness of the lower crust is variable and dependent on the inflow velocity of crustal materials. Here we build up both analytical and numerical models to study the mechanism and process of the uplift of Tibetan plateau at the eastern boundary. The results of the analytical model show that if the thickness of the lower crust can vary during the uplift process, the lower crustal viscosity of the Sichuan basin needs to be 1022 Pas to fit the observed elevation gradient. Such a viscosity is one-order magnitude larger than the previous results. Numerical model results further show that the state of stresses at the plateau boundary changes during uplift processes. Such a stress state change may cause the formation of different fault types in the Longmenshan fault area during its uplift history.  


Introduction
Tibetan plateau is the largest plateau in the world with an average elevation of *4500 m. The plateau has a flat interior and magnificent mountain belts around it (Tapponnier et al. 2001;Jiménez-Munt and Platt 2006) (Fig. 1). The mechanism for uplift of the Tibetan plateau has always been a focus in geological researches. The common consensus is that the uplift began when the Indian plate collided with the Eurasian plate after the oceanic lithosphere of the Tethys Sea subducted beneath the Eurasian plate (van der Voo et al. 1999;Royden et al. 2008). But the time of the collision is still under debate, probably starting from 50 Ma or earlier (Tapponnier et al. 2001;Yin and Harrison 2000;Royden et al. 2008;Wang et al. 2012). The suture zones distributed in Tibetan plateau imply that the plateau has undergone several episodes of different terranes being pieced together (Yin and Harrison 2000;Tapponnier et al. 2001). When being squeezed by the Indian plate, the crust of Tibetan plateau got shortened in south-north direction and thickened, with extensions at east-west direction and an outflow of materials toward east (Yin and Harrison 2000;Tapponnier et al. 2001;Zhang et al. 2004;Royden et al. 2008).
At the eastern margin of the plateau, the flow of materials breaks into two branches due to the obstruction of Sichuan basin: one toward Ordos basin and the other toward Indo-China block with a clockwise rotation Zhao et al. 2013). The elevation gradient varies significantly at the east boundary of Tibetan plateau. Close to the Sichuan basin, the elevation gradient is extremely large at the Longmenshan fault area, dropping from *4500 to *500 m within a short distance of *100 km. However, the elevation gradient becomes much gentler at the south and north sides of Sichuan basin D. Peng (&) Á W. Leng Laboratory of Seismology and Physics of Earth's Interior, School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China e-mail: pdd@mail.ustc.edu.cn ( Fig. 1). Understanding the uplift history and the elevation gradient variations at the eastern boundary of the Tibetan plateau is one of the key aspects for building up the dynamic evolution processes of the Tibetan plateau.
According to the results of seismic source mechanism, magnetotelluric imaging and anisotropic tomography results, the lower crust of the Tibetan plateau is probably a layer with high temperature and low viscosities (Nelson et al. 1996;Wei et al. 2001;Unsworth et al. 2005;Jiménez-Munt and Platt 2006;Priestley et al. 2008;Yang et al. 2012;Liu et al. 2014). First, the crustal thickness of the plateau is almost twice as much as the worldwide averaged value (Jiménez-Munt and Platt 2006). Therefore, it should contain more radioactive elements, leading to high temperature in the lower crust. Magnetotelluric results reveal that the electrical resistivity of the lower crust in the Tibetan plateau is lower than ambient regions such as Sichuan basin and Tarim basin, suggesting the possible existence of partial melting (Nelson et al. 1996;Wei et al. 2001;Unsworth et al. 2005). Second, according to the distribution of earthquake locations, the seismic activity is very weak in the lower crust (Priestley et al. 2008 Arrows with ellipses at the heads are the surface movement from GPS observation (Zhang et al. 2004), and the ellipses mean 95% error range. AA 0 notates the profile discussed later. The elevation of AA 0 is shown in Fig. 3. The inset map shows Tibetan plateau and ambient areas. The black lines are major plate boundaries. The blue box shows the location of our study area can be explained by the low mechanical strength of the lower crust. Tomography results also show low-velocity zones in the lower crust (Meissner and Mooney 1998;Galvé et al. 2006;Wang et al. 2007;Yao et al. 2008;Yang et al. 2012;Liu et al. 2014). Furthermore, small elastic thickness and high Q value of Tibetan plateau both support a weak lower crust for the Tibetan plateau (Jordan and Watts 2005;Wang et al. 2007;Bao et al. 2011;Zhao et al. 2013). Third, the lower crust of the plateau has strong seismic anisotropy, where the velocity of SH wave is much higher than SV wave (Shapiro et al. 2004;Huang et al. 2010;Yao et al. 2010). This radial anisotropy may be explained by the horizontal arrangement of anisotropic minerals such as mica (Shapiro et al. 2004;Huang et al. 2010;Shen et al. 2015). Such rearrangement of minerals may reflect flow motions in the lower crust layer and behaves as a direct evidence that the lower crust flow exists (Shapiro et al. 2004). Lower crustal channel flow may also be the reason for the extensional tectonics and, at the same time, the collapse of Songpan-Garze orogenic belt at early Mesozoic (Yan et al. 2008). Based on these observed factors of the lower crust, a crustal channel flow model has been proposed to explain the difference of elevation gradient at the east boundary of the Tibetan plateau (Royden et al. 1997;Clark and Royden 2000;Beaumont et al. 2001;Clark et al. 2005;Cook and Royden 2008;Royden et al. 2008). It proposed that the lower crust of the Tibetan plateau is mechanically weak with a lower viscosity than the upper crust and the upper mantle; therefore, the lower crust materials within the weak channel should be able to flow from the Tibet area to the ambient regions. In this model, the viscosity of the lower crust controls the elevations gradient at the plateau boundaries. A higher viscosity in the lower crust corresponds to a larger elevation gradient (Royden et al. 1997;Clark and Royden 2000). With this model and the observed elevation gradients, the viscosity of the lower crust beneath the Sichuan basin has been constrained as 10 21 Pas. However, when deriving this channel flow model, the thickness of the lower crust is kept constant. This is against the reality that the thickness of the lower crust should increase when new materials flow in. The thickness of the lower crust plays an important role in surface topography and elevation gradient. It is essential to build a more realistic crustal flow model to incorporate this effect.
Here we study the uplift process at the Tibet-Sichuan boundary, so we choose a profile across the Longmenshan fault (profile AA 0 in Fig. 1). In this study, we first analytically derive an improved channel flow model in which the thickness of the lower crust is controlled by the inflow velocity of the crustal materials. Then with the observed elevation gradients, this new model is used to constrain the viscosity of the lower crust beneath the Sichuan basin. Starting from this analytic model, we further build numerical models to explore the time evolutions of the uplift history in the Longmenshan fault area. We also discuss the temporal and spatial distributions of stresses in this area during the uplift processes.

An improved analytical channel flow model
During the strong uplift processes of the Tibetan plateau, no significant crustal shortening was observed (Yin and Harrison 2000;Clark and Royden 2000). Therefore, Clark and Royden (2000) proposed an analytical channel flow model to explain the uplift mechanism, in which weak lower crust materials can horizontally flow within a channel of constant thickness (Fig. 2a). When the flowing materials are blocked by strong basins such as the Sichuan basin, lower crustal materials accumulate and cause increased crustal thickness which leads to high elevations. The elevation gradient is dominantly controlled by the lower crust viscosity of the strong basins. With the observed elevation gradient at the Tibet-Sichuan boundary, close to the Longmenshan area, Clark and Royden (2000) constrain the viscosity of the lower crust beneath the Sichuan basin to be 10 21 Pas. There are two points in the previous channel flow model which need to be reconsidered. First, the thickness of the channel is fixed as a constant (Fig. 2a), which is obviously not in accordance with the tectonic facts that the thickness of the flow channel should gradually increase with the influx materials (Fig. 2b). Second, Clark and Royden (2000) chose 20 Ma as the duration time for the uplift process at the Longmenshan area. But Kirby et al. (2002) used thermochronology to suggest that the uplift of Longmenshan started at 11-5 Ma. Wang et al. (2012) also proposed that there may be two rapid uplift processes at the Longmenshan area, and the later one probably started from 15 to 10 Ma.
Here we set up an improved channel flow model to incorporate the variations of the channel thickness and the updated duration time for the uplift process at the Longmenshan area. In our model, the uplift of Longmenshan area began at 11 Ma. The channel thickness is a variable which depends on the amount of materials flowing in from the left boundary (Fig. 2b). The velocities above and below the channel are zero. The left and right boundaries are free slip. The thickness of the lower crust is set as 27 and 15 km at these two boundaries, which, respectively, represent the lower crust thickness beneath the Tibetan plateau and the Sichuan basin from the P-wave velocity studies (Wang et al. 2007) (Fig. 2b). The initial thickness of the channel is 15 km. The crust and mantle density is 2600 and 3300 kg/m 3 , respectively. The viscosity of the lower crust materials in the channel is taken as a constant.
In this model, although the upper and lower boundaries of the flow are not parallel, we can still consider flow in the channel as a Poiseuille flow in each infinitely small box in the horizontal direction. According to the Poiseuille law, the horizontal velocity u depends on the horizontal pressure gradient (Turcotte and Schubert 2002): where x and z are the horizontal and vertical coordinates; l is viscosity. p ¼ q c gE is pressure, where q c , g and E are the crustal density, gravitational acceleration and surface topography. h is the thickness of the lower crust channel. The total material flux in the channel, U, is: The temporal variation of lower crust thickness c is the same as the total material flux: The surface elevation is controlled by crustal isostasy as If we consider the thickness of the channel h is the same as the thickness of the lower crust c, we finally obtain the controlling equations for channel thickness and the surface topography: where a ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; q m is the mantle density. Equations (5) and (6) are solved with Wolfram Mathematica to obtain h and E at any instant time when initial and boundary conditions are given. According to the crustal structure from Wang et al. (2007) and topography data, We plot the temporal variations of lower crustal thickness and surface elevation at x = 10, 50 and 100 km (Fig. 3a). It can be observed that near the left boundary (red lines in Fig. 3a), both the thickness of the lower crust and the elevation increase quickly in the first few million years. Then they gradually  Fig. 3a), the thickness and elevation grow much slower. As pointed out by Clark and Royden (2000), the elevation gradient at the Tibet-Sichuan boundary is controlled by the viscosity of the lower crust beneath the Sichuan basin. The larger the lower crust viscosity, the steeper the elevation gradient. In order to reproduce the observed elevation gradient, the lower crust viscosity is constrained to be 10 22 Pas (Fig. 3b). This value is one-order magnitude larger than previous modeling results, suggesting that a variable channel thickness plays an important role in the channel flow model. When the lower crust channel is thickened, more materials can flow into Sichuan basin through the channel. Therefore, to form the same elevation gradient, we need a larger viscosity of lower crust in our model.

Numerical model results
The above analytical model provides us insights into the channel flow evolution, and the uplift processes with a variable channel thickness provide a more realistic lower crust viscosity 10 22 Pas. But it is limited for parameter choices and tectonic setting due to model simplicity. For example, the upper crust is ignored, and the material is purely viscous fluid. Particularly, we are interested in the stress states in the upper crust during the uplift processes, which can only be obtained from numerical modeling. Therefore, we further develop a numerical model to study the uplift process at the Tibet-Sichuan boundary. We use a 2D viscoelastic model to study the crustal thickening and uplift processes. The governing equations for the model can be referred to our previous papers Gurnis 2011, 2015). Here we only introduce model geometry and parameters.
We model the channel flow with a 360 km by 120 km rectangular box in which the temperature effects are ignored. The model is divided into four layers from top to the bottom: sticky air layer, upper crust, lower crust and the mantle (Fig. 4). The sticky air layer is compressible and used for simulating a free surface (the same as in Leng and Gurnis 2015). The thickness and density of the upper crust, lower crust and lithospheric mantle are set as 28, 15 and 62 km, and 2730, 2980 and 3300 kg/m 3 , respectively. The Sichuan basin lower crust viscosity is taken from the analytical model, which is 10 22 Pas. After a series of model tests, the viscosities for the upper crust and the lithospheric mantle are fixed at 10 24 and 10 23 Pas. In the horizontal direction, we consider three different types of models according to the viscosities in the lower crust channel: uniform model (Case 10, Fig. 4a), sharp-edge model (Case 20, Fig. 4b) and transitional edge model (Case 30, Fig. 4c).
The uniform model represents the Sichuan basin with a lower crust viscosity of 10 22 Pas. The sharp-edge model is divided into two parts at x = 120 km. The right part represents the Sichuan basin, and the left part represents the Tibetan plateau with a smaller lower crust viscosity of 10 21 Pas. For the transitional edge model, an extra block is inserted between x = 60 and 120 km. This block has a lower crust viscosity of 5 9 10 20 Pas, which represents a weak Longmenshan fault area as suggested by previous studies (Zhang et al. 2009). All the model boundary conditions are free slip, except for the left boundary of lower crust where an inflow boundary condition is used. The model resolution is 129 9 65 grids.
In the analytical model, the channel flow is driven by the horizontal pressure gradient. Here in the numerical model, we add an influx boundary condition at the left boundary to simulate the material transport. Then the magnitude of the influx velocity needs to be determined. For the uniform model (Fig. 4a), we test three different influx velocities, 0.39, 0.66 and 1.31 cm/a. Model results show that after *11 million years of uplift process, the surface elevations at the left boundary are 3000, 4000 and 7000 m, respectively. We, therefore, choose 0.66 cm/a, which causes a surface elevation consistent with the Tibetan plateau, as the reference influx velocity for all our models. Such a velocity also agrees to the velocity difference between the Tibetan plateau and the Sichuan basin from GPS data (Zhang et al. 2004).
With the influx velocities at the left boundary, the lower crust is thickened with time. Figure 5 shows the viscosities, velocities, strain rates and principal stresses for Case 10 after 11 million years of evolution. It can be observed that velocities decrease quickly from the model boundaries to the interior (Fig. 5b). The normal strain rates in the horizontal and vertical directions (e 11 and e 33 ) both reach the maximum values in the lower crust close to the left boundary ( Fig. 5b and 5c), indicating that the deformation is most intense in this region. But due to the low viscosity in the lower crust, the principal stresses in the lower crust are generally smaller than those in the upper crust and the lithospheric mantle (Fig. 5c). Figures 6 and 7 show the viscosities, velocities, strain rates and principal stresses for Case 20 and Case 30, respectively. The thickness of lower crust in Case 20 (Fig. 6a) and Case 30 (Fig. 7a) is smaller compared with Case 10 (Fig. 5a), which indicates lower topography, while the strain rate in Case 20 (Fig. 6b, c) and Case 30 (Fig. 7a, c) is larger than Case 10 (Fig. 5b, c). Because of the smaller viscosity at the left boundary, the stress in Case 20 (Fig. 6c) and Case 30 (Fig. 7c) is much smaller than that in Case 10 (Fig. 5c).
We first look at the uplift processes for the uniform model Case 10. The temporal evolution of surface elevations at three locations (x = 10, 50 and 100 km) is plotted versus time (solid lines in Fig. 8). It can be observed that the elevation at x = 10 km increases with time and reaches 4000 m after *10 Ma. The elevations at x = 50 and 100 km increase much more slowly with time and therefore form a large elevation gradient in this uplifted area. For the sharp-edge model Case 20 (dashed lines in Fig. 8), elevation increases more slowly at x = 10 km but faster at The magnitude and orientation of surface stress determine the fracture formation and fault types in the uplifted zone. We then investigate the surface stress states for these three different models. For Case 10, the magnitude of principal stress decreases with time at x = 10 and 100 km, but keeps almost constant at x = 50 km (solid lines in Fig. 9). The orientation of the principal stress can be expressed with angles between the extensional axis of the stress and the x direction. For Case 10, the compressional axis at x = 10 and 100 km keeps to be horizontal and vertical for the whole uplift process, respectively (red and green solid lines in Fig. 10). This suggests that the crust undergoes strong horizontal extension in high-elevation regions but strong compression in low-elevation regions. At x = 50 km, the extensional axis gradually rotates toward horizontal direction during the uplift process (blue solid line in Fig. 10). For Case 20 with a viscosity jump in the lower crust, the magnitudes of principal stress also decrease with time at x = 10 km and keep almost constant at x = 50 km (red and blue dashed lines in Fig. 9). But at the low-elevation region where x = 100 km, it can be observed that the magnitudes of stress first decrease to be     almost zero at *9 Ma and then start to increase again (green dashed lines in Fig. 9). More interestingly, the orientations of the extensional axes rotate by almost 90 degrees at the same time (green dashed lines in Fig. 10). This suggests that the stress state at this location undergoes an orientation transition from strong compression to strong extension in the horizontal direction. Such a transition in stress state may cause different fault types in the Longmenshan area during the uplift processes. For Case 30, the surface stress states (dotted lines in Fig. 9) are similar to the results of Case 20, but the orientations of the extensional axes behave differently (dotted lines in Fig. 10). For Case 30, the rotation appears *2 Ma ahead of that in Case 20, which suggests that if the extra weak block exists between Tibetan plateau and Sichuan basin, the transition of fault types will happen earlier.

Discussion
Except for the channel flow model, there are other models to explain the uplift of the Tibetan plateau. Tapponnier et al. (1982Tapponnier et al. ( , 2001 considered that the Tibetan plateau contains several rigid blocks, which were broken down and rotated because of the subduction of Indian plate. But GPS data of surface motion suggest that the deformation process is more continuous than rigid (Zhang et al. 2004;Yang et al. 2012). Molnar et al. (1993) thought that at about 10 Ma, a rapid uplift happened which raised the elevation by more than 1000 m. To explain this, they came up with a gravitational instability theory. It proposed that after the crust is thickened, lithosphere root goes deep into the asthenosphere. Then mantle convection removes the bottom of lithosphere and takes its place with hot lighter asthenosphere material, leading to surface uplift. Jiménez-Munt and Platt (2006) agreed to this opinion and ran several numerical models based on it. Their model can also explain the flat interior of Tibetan plateau and the extension in east-west direction. It is still in debate which one, the channel flow model or the gravitational instability model, is the real mechanism for the uplift of Tibetan plateau.
We used a 2D model to study the uplift processes at the Tibet-Sichuan boundaries. But at the east boundary of the plateau, the material flow breaks into different branches toward Ordos basin and Indo-China block. A full 3D model is needed to better reveal the effects of the ambient regions and different elevation gradients along the whole east boundaries of the plateau. On the other hand, the S wave anisotropy is revealed by several studies in this region (Shapiro et al. 2004;Huang et al. 2010. Such kind of S wave anisotropy also needs a 3D numerical model to fully simulate the corresponding flow pattern. Additionally, when elevation gradient is large, the effect of erosion may be very significant. Beaumont et al. (2001) have designed models for the Himalayan tectonics which combine lower crustal channel flow and surface denudation. The outcrops of lower may indicate the extension and thinning of the crust and thus reduce the topography (Yan et al. 2008). Both these processes can retard the uplift of topography, and their effects need to be further studied.
We did not consider the effects of large faults on the uplift processes. Active faults are widely distributed along the plateau boundaries, among which there are large ones such as Longmenshan fault, Red River fault, Xianshuihe fault. These large fault zones may also play important roles in the uplift processes of the plateau boundary. The effects of these faults should be considered and analyzed in future models.

Conclusions
In summary, we build an improved analytical channel flow model in which the thickness of the lower crust channel can vary with time according to the amount of materials flowing into the channel. With such a model and observed elevation gradient at the boundaries between the Tibetan plateau and the Sichuan basin, the lower crust viscosity of the Sichuan basin is constrained as 10 22 Pas, which is one-order magnitude larger than the previous model with a constant channel thickness. Moreover, we use this lower crust viscosity and build three different types of numerical models to study the uplift process and stress state at the boundaries between the Tibetan plateau and the Sichuan basin. Our results show that when there is a viscosity jump in the lower crust channel, the stresses at the surface in the Longmenshan fault area could rotate from horizontal compression to extension during the uplift process. Such a stress state transition may lead to the formation of different fault types in the uplifting process of the Longmenshan area.