Influence of Regional Erosion and Sedimentary Loading on Fault Activities in Active Fold-Thrust Belts: Insights From Discrete Element Simulation and the Southern and Central Longmen Shan Fold-Thrust Belt

Four groups of discrete element models (DEMs) were set-up to simulate and analyze the influence of regional erosion and sedimentary loading on the formation and spatial-temporal evolution of faults in the southern and central Longmen Shan (LMS) active fold-thrust belt. The interior characteristics of faults in the southern and central LMS fold-thrust belt were also evaluated during the interaction of tectonic processes and surface processes according to the stress-strain analysis from DEM results. The results showed that synkinematic erosion promoted the reactivation of pre-existing faults in thrust wedges and also retarded the formation and development of new incipient faults in the pre-wedge regions. Meanwhile, synkinematic sedimentation also delayed the development of new incipient faults in the pre-wedge regions by promoting the development of thrust faults in the front of thrust wedges, causing these thrust wedges in supercritical stages with relatively narrow wedge lengths. According to these DEM results, we infer that: 1) The characteristics of erosion and sedimentation in the central and southern LMS have important influences on the activities of large faults which are extended into the deep detachment layer; 2) Besides differential erosion, the differential sedimentary loading may also be one of the important factors for the along-strike differential evolution of the LMS fold-thrust belt. This kind of differential deposition may lead to differential fault activity and uplift in the interior thrust wedge and pre-wedge region in the central and southern LMS; 3) Compared to the northern LMS, the central LMS and southern LMS is more conducive to the occurrence of earthquakes, because of synkinematic sedimentation (such as the growth of Chengdu plain) has a greater blocking effect on the stress propagation and strain convergence on the fault planes of front faults of an active thrust wedge.


INTRODUCTION
The Longmen Shan (LMS) Mountains, located in the eastern margin of the Qinghai-Tibetan Plateau with an NE trend, was a suture zone between the Yangtze block and the Songpan-Ganzi block and was reactivated due to the far-field compression of the India-Asia collision (Burchfiel et al., 1995;Roger et al., 2003, ;Zhou et al., 2006;Roger et., 2010;Yan et al., 2011) (Figure 1A). It has experienced a long history of deformation since the Mesozoic, which resulted in the development of a series of fold-thrust belts along the western boundary of the Sichuan Basin (Burchfiel et al., 1995;Yan et al., 2008), and the deposition of a suite of late Triassic to Cretaceous-Quaternary sediments in the Sichuan Basin Li Z. et al., 2018;Luo et al., 2021) (Figures 1B,C). The 2008 M W 7.9 Wenchuan earthquake occurred deep in the interior of the LMS fold-thrust belt, 200 km northwest of the toe of the fold-thrust belt in the Sichuan basin , accompanied by the formation of coseismic reverse-and oblique-slip surface faulting by two large thrust faults (Lin et al., 2009;Liu-Zeng et al., 2009;Xu et al., 2009;Li et al., 2010). The occurrence of this earthquake showed that the LMS fold-thrust belt is still active nowadays and is characterized by strong tectonic uplift, rapid erosion, and devastating seismic hazards (Hubbard and Shaw, 2009;Xu et al., 2009;Yin, 2010;Gao et al., 2016;Tan et al., 2017;Wang et al., 2021).
For many years, many scholars have been engaged in the research on the tectonic geomorphology, tectonic activities and surface processes of LMS with various methods. The LMS foldthrust belt has already been a natural laboratory for studying the correlation between surface processes (such as erosion system, sedimentary system and climate system) and tectonic processes (such as fault activity and tectonic uplift). Especially in recent years, the interaction between surface and tectonic processes in the LMS fold-thrust belt has become a research hotspot. For example, Gao et al. (2016) applied two geomorphic indices to evaluate the differential uplift of the central and southern LMS fold-thrust belt, which suggested the relationship between the regional erosion and the rates of tectonic uplift. Tan et al. (2018) indicated that external forces play a dominant role for the alongstrike variation of fault activity and in thrust belt evolution, such as the fluvial incision (erosion), which may further produce the different rupture behavior of the coseismic slip partitioning onto two sub-parallel faults along the LMS fold-thrust belt during the 2008 M W 7.9 earthquake. Liu et al. (2020) highlighted the role of erosion in creating thrust recesses in a critical-taper wedge based on mechanical models. Wang et al. (2021) emphasized the role of tectonic structures in regulating the pattern of denudation and topography across the LMS fold-thrust belt. All these studies mentioned above indicated that denudation played an important role on the structural evolution of the LMS fold-thrust belt. However, the effects of sedimentary loading were less discussed. More and more studies have shown that sedimentary loading also played important roles on the evolution of orogenic fold-thrust belts, e.g., DeCelles and Carrapa (2021) have studied the deformation-erosionsedimentation interactions in a salient-recess pair in the Himalayan thrust belt, and Liu et al. (2021) have shown effects of deposition on fault activity (e.g., the effects on constraining vertical slip rates of thrust faults) in the north Qilian Shan, NE Tibetan Plateau. Moreover, in the LMS foldthrust belt, Li et al. (2016), Li C. et al. (2018) identified growth strata and revealed early to late Pleistocene activity on the Range Front blind thrust. Recently, by using sandbox analog experiments, Luo et al. (2021) investigated the influences of differential rates of synkinematic erosion and sedimentation on wedge geometries and fault activities during the development of thrust wedges, and further discussed the influences from the change of erosion-sedimentation along the strike on the structural evolution of the LMS fold-thrust belt; Mao et al. (2021) have shown that surface erosion and sedimentation affected the position, morphology, lateral propagation, and connection of new structures in fold-thrust belts. Overall, sufficient studies have shown the important relationship between surface processes (denudation and/or deposition) and tectonic processes (fault activity) in fold-thrust belts, although there are still some different viewpoints on the main controlling factors. Therefore, how do surface processes, such as synkinematic erosion and sedimentary loading, influence the development of a fold-thrust wedge, in particular the characteristics of stress and strain in the interior wedge of the LMS fold-thrust belt? It is worthy of further research from the perspective of quantitative simulation experiments.
In this paper, based on previous studies, we simulate, discuss, and analyze the influences of synkinematic erosion and sedimentary loading on the formation and spatial-temporal evolution of faults in the southern and central LMS active fold-thrust belt from the perspective of numerical modeling (discrete element model, DEM) experiments. The interior characteristics of faults in the southern and central LMS foldthrust belt are evaluated during the interaction of tectonic processes and surface processes according to the stress-strain analysis from DEM results. In addition, our modeling results have the experimental reference significance on the study of the dynamic background of tectonic uplift in the LMS foldthrust belt.

GEOLOGICAL SETTING
The LMS is an active fold-thrust belt, composed of a series of NE trending thrust faults between the Tibetan Plateau and the Sichuan basin ( Figure 1B) (Avouac and Tapponnier, 1993;Xu and Kamp, 2000;Wang et al., 2012;Li et al., 2016). The lithology and fault characteristics of the LMS fold-thrust belt were determined by means of field geological investigation and remote sensing interpretation (e.g., Deng et al., 1994;Shen et al., 2009;Zhang et al., 2011;Gao et al., 2013;Ren et al., 2014 and many others), geophysical data (such as industrial seismic profiles) and drilling data (e.g., Chen et al., 2005;Jia et al., 2006;Hubbard and Shaw, 2009;Hubbard et al., 2010;Li et al., 2010;Li Z. et al., 2013;Lu et al., 2014;Wang et al., 2014;Li et al., 2016;Li C. et al., 2018). Especially after the 2008 Wenchuan earthquake, a lot of field investigation works have been carried out for surface ruptures (e.g., Wang et al., 2008;Hubbard and Shaw, 2009;Lin et al., 2009;Liu-Zeng et al., 2009;Shen et al., 2009;Xu et al., 2009;An et al., 2010;Feng et al., 2010;Hashimoto et al., 2010;Li et al., 2010;Lin et al., 2010;Yin, 2010). Four major imbricate thrust systems were interpreted consistently : The Back range fault system, Central range fault system, Front range fault system, and Front blind range fault system. Correspondingly, several thrust faults can be interpreted from NW to SE in these imbricate thrust systems. There are the Wenchuan-Maoxian fault (WMF) in the Back range fault system, the Yingxiu-Beichuan fault (YBF) and Wulong Fault (WLF) in the Central range fault system, the Pengxian-Guanxian fault (PGF) and Shuangshi-Dachuan Fault (SDF) in the Front range fault system, and the range front blind thrust (RFBT) in the Front blind range fault system (Figures 2, 3).
For the tectonic analysis of the southern and central LMS, due to the lack of subsurface controls on the fault geometry, there are some differences in the interpretation of deep structures. Hubbard and Shaw (2009) inferred that the main detachment is close to the depth of 20 km according to the structural analysis of the balanced geologic cross-section (Figure 2), whereas another research group Li et al., 2010;Li et al., 2016) suggested that the YBF and PGF sole into the main detachment at a depth of 15-17 km based on seismicreflection profiles, drill-hole data, and relocated aftershocks of the 2008 Wenchuan earthquake (Figure 3). In addition, according to the results of structural analysis and modeling, since the late Cenozoic, the total shortening displacement of the southern and central LMS fold-thrust belt can reach 20-40 km, with a shortening ratio of 19.3-39.2% (Chen et al., 2005;Hubbard and Shaw, 2009;Hubbard et al., 2010;Li Z. et al., 2013). Moreover, according to Hubbard et al. (2010), further information about the LMS fold-thrust belt was obtained: 1) Upper-crustal shortening plays an important role in developing the topography of the LMS fold-thrust belt; 2) Two taper wedges (two differential topographic slopes in the front of the LMS and the Sichuan basin, respectively) could be defined from the hinterland to foreland: A steep surface slope caused by a slightly strong basal detachment in the front of the LMS, and a shallow surface slope caused by a middle level, weak detachment in the Sichuan basin; 3) A minimum of 95% shortening was absorbed in the frontal LMS fold-thrust belt according to the retrodeformed cross sections.
In recent years, many important achievements have also been published in this study of erosion rates along-strike or cross-strike in major faults (Xu and Kamp, 2000;Tan et al., 2015;Gao et al., 2016;Tan et al., 2017;Tan et al., 2019;Luo et al., 2021;Wang et al., 2021), as well as the syntectonic sedimentation deposited in the LMS active fold-thrust belt (Yan et al., 2013;Li et al., 2014;Li et al., 2016;Li C. et al., 2018;Luo et al., 2021). According to previous research on the regional erosion system, at least 4-5 km of strata have been eroded in the southern and central LMS at an erosion rate of 0 mm/yr to 1.0 mm/yr from the range front to the hinterland of the LMS foldthrust belt since the late Miocene (Tan et al., 2017;Shen et al., 2019;Ye, 2021;Wang et al., 2021). In addition, the differential erosion characteristics and syntectonic sedimentary loading along-strike of the LMS fold-thrust belt were further confirmed. According to the analysis of erosion rates by thermochronology data (e.g., Richardson et al., 2008;Li et al., 2012;Cook et al., 2013;Richardson et ., 2013;Tan et al., 2015, Tan et al., 2017Jia et al., 2020;Ye, 2021) and the thickness of syntectonic sedimentary strata by field exploration data (such as drilling data) Li Z. et al., 2018;Chen et al., 2019), it is presented that the erosion rates and the thickness of syntectonic sedimentary strata in the central and southern LMS fold-thrust belt are heterogeneous. Taking the fault of PGF or SDF as the dividing line in the piedmont fault zone, the erosion rates of the footwall and hanging wall strata are respectively about 0.3-0.4 mm/yr and 0.5-0.6 mm/yr since ∼6-8 Ma in the southern LMS ( Figure 2); while in the central LMS, the erosion rates of the footwall and hanging wall strata are respectively about less than 0.1 mm/yr and 0.2-0.5 mm/yr since ∼7.5 Ma ( Figure 3).
During the Late Triassic to Early Jurassic, the Songpan Ganzi fold belt was thrust upon the western margin of the Yangtze Plate (Yan et al., 2011). The LMS fold-thrust belt uplifted rapidly and was accompanied by approximately 1-2 km of wedge-shaped molasse deposited in the Chengdu plain, an alluvial plain in the western Sichuan basin ( Figure 1B) (Luo et al., 2021). Moreover, the Cenozoic sedimentary strata, including the Eocene-Oligocene Minshan and Lushan formations, the upper Pliocene Dayi Formation, the middle Pleistocene Ya'an Formation, the upper Pleistocene Chengdu Clay Formation, and uppermost Holocene deposits, crop out mainly on the Chengdu plain and can also be found locally in the foothills of the LMS . The thickness of Cenozoic syntectonic sedimentary strata in the piedmont of the southern and central LMS (southern Chengdu plain and northern Chengdu plain, respectively) is about 0-550 m ( Figure 1C). The depocenter is located in front of the central LMS, which caused the maximum sedimentary thickness in the central LMS to be almost 200 m thicker than that in the southern LMS ( Figure 1C). These strata have recorded the interactions between tectonic deformation and sedimentation along the LMS range front ( Figure 2B and Figure 3B).  Figure 1; data of average erosion rate (AER) reference from Ye, 2021; epicenter data from Li et al., 2010); (B) a geological profile shows the synkinematic sedimentation deposited in the southern LMS active fold-thrust belt (Li Z. et al., 2018). The profiles show that there are typical characteristics of synkinematic erosion and sedimentary loading in the southern LMS fold-thrust belt.

MODEL SET-UP
The Discrete Element Modeling (DEM), derived from the molecular dynamics method where the interaction forces between particles are set ( Figure 4), allows large displacement inside the models (Wu et al., 2019). In this model, strata are characterized as independent elastic particles of different sizes. It is suitable to simulate dynamic fault formation and evolution by  Figure 1; data of average erosion rate (AER) reference from Tan et al., 2017). The measured coseismic slip of the surface rupture  indicates that the displacement in the central segment is dominated by thrust slip. (B) Interpretation of the shallow seismic profile that exhibits the presence of the synkinematic strata. These profiles also show that there are typical characteristics of regional erosion and sedimentary loading in the central LMS fold-thrust belt.
FIGURE 4 | Schematic diagram of particle interactions in DEM. (A) the interaction between the boundary wall and particles of strata, which cause particles to move, rotating by overlapping with each other under the role of self-gravity and external stress transmitted from the moving boundary wall; (B) two forces (normal force and shear force) during the interaction between two particles.
In this paper, four two-dimensional DEM models (Table 1) are designed with an initial set up of 60 units' length and 6 units' height. 24,309 particles with two different diameters are generated randomly ( Figure 5). The parameters of particles are shown in Table 2. Two different types of rock particles are defined by setting interparticle bond properties among particles (shown in Table 3), where pre-shortening strata (brittle strata), syn-shortening sedimentation (brittle strata), and detachment stratum (ductile strata) are shown as black and yellow layers, middle gray layer, and red layer, respectively. The cohesion and the friction angle of brittle strata are ∼10.5 MPa and ∼18.6, respectively, which were obtained through two-dimensional biaxial compression (Morgan, 2015). In addition, plate walls are shown in purple ( Figure 5). The particles' friction coefficient of brittle strata and the plate walls is 0.3. According to the previous research on the friction coefficient of basal detachment in the LMS thrust wedge , the slightly strong value ( 0.1) of the friction coefficient of basal detachment is selected in our models. Moreover, there is no cohesion between the particles of the detachment and the plate walls in our models (Morgan, 2015).
These four DEM experiments were run to analyze the influences of erosion and sedimentary loading on the activity of the range front faults in the southern and central LMS foldthrust belt. In our models, 1 unit represented 2 km. The total shortening displacement was 10 units, which represented the regional shortening displacement of 20 km in nature, or a   (Li et al., 2017;Li C. et al., 2018;Li, 2019; was used in our experiments. To analyze the distribution characteristics of stress and strain in the fault zone in our three DEM experiments, stress and strain were calculated based on the method proposed by Morgan (2015).  Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 659682 7

MODEL RESULTS
In this study, four representative experiments ( Table 1) were analyzed. The baseline model (Model 1) had no synkinematic erosion or sedimentation. Models two and three only had synkinematic erosion and sedimentation, respectively. Model four had both synkinematic erosion and sedimentation. The structural evolution of these four DEM models is described in detail and shown in

Structural Evolution
The experiment results show that a back thrusting fault occurred near the backstop according to the boundary effect of models in all models. Moreover, differential evolution characteristics of their thrusting faults are presented (shown in Figures 6, 7) among these four DEM models by altering experimental  Figure 6A shows the progressive evolution of the baseline Model one in six stages up to a maximum of 10 units of shortening. In this model, the initial 2 units of shortening (Panel a2 in Figure 6A) formed a forward-breaking thrust fault F1 together with an incipient back thrust and its associated flat-topped hanging-wall anticline. A low-relief detachment fold formed a gentle bulge and caused local thickening of the shallow stratain front of the thrust fault F1, marking the position of the incipient thrust fault F2 (Panel a2 in Figure 6A). By 4 units of contraction thrust, F2 had propagated to the surface forming a frontal ramp fold (Panel a3 in Figure 6A).

Model 1 With no Synkinematic Erosion or Sedimentation
Subsequent shortening shown in Panels a4-a6 in Figure 6A produced an imbricate thrust system including a total of four forward-vergent thrust faults F1-F4, developed in turn with differential fault slip displacements ( Figure 7A). During this shortening stage, when the next thrust fault starts to develop, the previous one basically stops developing. The model shortening is mainly absorbed by the thrust faults F2 and F4 ( Figure 7A). After a total of 10 units of shortening, the final model consisted of a thrust wedge with a forward-sloping surface topography formed by four foreland-vergent thrust faults each with their associated hanging-wall ramp folds (Panel a6 in Figure 6A). These thrust faults had nucleated in a forwardbreaking sequence. The location of the front surface slope toe of this thrust wedge in model one was at 30 units from the backstop ( Figure 6A) with a taper angle of 7.2°( Figure 8A).

Model 2 With Synkinematic Erosion and no Sedimentation
In model 2 ( Figure 6B), the initial geometries and properties of the model evolution (Panels b1-b3 in Figure 6B) were identical to those in Model one described above (Panels a1-a3 in Figure 6A). Then different from Model 1, subsequent shortening increments, beyond the initial 5 units of contraction, were accompanied by the introduction of synkinematic erosion by removing the strata above 6.5 units from the top surface of the thrust wedge composed by two forward-breaking thrust faults F1 and F2 (Panel b4 in Figure 6B).
With continued shortening after this erosion, the thrust fault F2 continued its activity ( Figure 7B), while in Model one it was almost stopped to develop at this stage because of the formation of new thrust fault F3 ( Figure 7A). Until 6 units of shortening, the thrust fault F3 began to develop strongly in this model (Panel b4 Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 659682 9 in Figure 6B and Figure 7B). Subsequent shortening shown in Panels b5 and b6 in Figure 6B indicated that thrust fault F3 still remained strongly active until the occurrence of thrust fault F4 ( Figure 7B). By a shortening of 10 units, the slip displacements of faults F2 and F3 are significantly larger than those in Model 1, while the slip displacement of fault F4 in Model two is significantly smaller than that in Model 1 (Figures 7A,B). The development of pre-existing faults (F2 and F3) in Model two was promoted after erosion, and the development of new thrust fault F4 was inhibited. Therefore, the location of the surface slope toe of this thrust wedge in Model two was at 28 units from the backstop ( Figure 8B) with a taper angle of 8.8°( Figure 8A), which is shorter than that in Model 1.

Model 3 With Synkinematic Sedimentation and no Erosion
In Model 3 ( Figure 6C), the initial geometries and properties of the model evolution ( Figures 6C1-C3) were also identical with those in Model one described above (Panels a1-a3 in Figure 6A). After 5 units of shortening, a syn-kinematic sedimentary layer was considered based on Model one in the front of the thrust wedge composed of two forward-breaking thrust faults F1 and F2 (Panel c3 in Figure 6C). After 6 units of shortening, a new forward-breaking thrust fault F3 formed with a larger fault slip displacement than that in Model 1 (Panel a4 in Figure 6A, Panel c4 in Figure 6C, and Figures 7A,C). However, during this stage, compared to the incipient development of the thrust fault F4 in Model 1 (Panel a4 in Figure 6A), no features indicated that a new thrust fault was about to develop in Model 3 (Panel c4 in Figure 6C). Up to 10 units of shortening, the forwardbreaking thrust fault F4 began to develop in Model 3. The development of this new thrust fault was inhibited obviously by sedimentary loading compared with Model 1 (Panel a6 in Figure 6A and Panel c6 in Figure 6C). The development of new thrust faults in Model three was delayed comparing with the previous two models, which caused the location of the surface slope toe of thrust wedge in Model three to be only at 24 units from the backstop ( Figure 8B) with a taper angle of 13.9°( Figure 8A). The model shortening was mainly absorbed by the thrust faults F2 and F3 ( Figure 7C).

Model 4 With Synkinematic Erosion and Sedimentation
In Model 4 ( Figure 6D), the initial geometries and properties of the model evolution (Panels d1-d3 in Figure 6D) were also identical with those in Model one described above (Panels a1-a3 in Figure 6A). After an initial shortening of 5 units, the hinterland of the wedge was eroded by removing the strata above 6.5 units all at once from the top surface of the thrust wedge composed by two forward-breaking thrust faults F1 and F2 (Panel d4 in Figure 6D). A synkinematic layer was also deposited onto the front slope of the growing thrust wedge after erosion (Panel Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 659682 d4 in Figure 6D). After 6 units of shortening, the thrust fault F2 still kept on activity ( Figure 7D), and the thrust fault F3 began to develop in this model (Panel d4 in Figure 6D). Until 7 units of shortening, the thrust fault F3 began to propagate strongly ( Figure 7D). After 10 units of shortening, the thrust fault F3 was still keeping active and cutting through the sedimentation layer (Panel d6 in Figure 6D). There was no new thrust fault occurring in front of the thrust fault F3 (Panel d6 in Figure 6D and Figure 7D). The development of thrust faults F1-F3 in this model caused the location of the surface slope toe of the thrust wedge in Model four to be at 23 units from the backstop ( Figure 8B) with a taper angle of 10.9°( Figure 8A). The model shortening was mainly absorbed by the thrust faults F2 and F3 ( Figure 7D).
Critical Wedge Geometries Figure 8A shows the incremental wedge geometries as indicated by their taper angle, α + β, plotted against horizontal shortening. In these four models, the taper angle was measured as the frontal topographic slope α, because the dip of the basal detachment (β) was set to 0°. According to the previous research on the taper angle in a fold-thrust belt from Davis et al. (1983), McClay (2011), andYang et al. (2017) by theoretical calculation, analogue modeling experiments, or numerical modeling experiments, the critical taper should be 9°-10°, which is almost consistent with our DEM models in this paper ( Figure 8A). In our models, the wedge geometries had high slopes (greater than 20°) for the first 5 units of shortening as the backstop geometries were formed ( Figure 8A). Beyond the first 5 units of shortening, Model 1 with no synkinematic erosion or sedimentation presented a subcritical stage along with the development of thrust faults F3 and F4 after 7 units of shortening ( Figure 8A). Model 2 with synkinematic erosion displayed a cyclic behavior with taper angles oscillating above and below the critical taper angle ( Figure 8A). In Model 3 with synkinematic sedimentation and no erosion and Model 4 with synkinematic erosion and sedimentation, the thrust fault F3 in the thrust wedge, remaining always active after 5 units of shortening ( Figures 7C,D), meant the thrust wedge was kept in the supercritical stage, respectively ( Figure 8A).
The cyclic nature of the thrust wedge geometries and the distribution of the internal deformation within the overall thrust wedge is also reflected in the plot of wedge length as indicated by the distance between the deformation front in the foreland and the initial location of the backstop ( Figure 8B). Significant increases in the lengths of the imbricate thrust wedges occurred each time a new forward-breaking thrust fault developed. Propagation of the deformation front into the foreland was inhibited by synkinematic sedimentation and is reflected by an increased slope in the curves for Model three and 4 ( Figure 8B). Compared with Model one and Model 2, the interval time of the episodic motion on particular high-displacement thrusts was extended in Model three and Model 4.
We also plotted the wedge height against the horizontal shortening ( Figure 8C). Comparing Model three to Model 1, we found that the wedge height of Model 3 with synkinematic sedimentation is higher than that of model one without synkinematic sedimentation. Moreover, compared Model four to Model 2, the wedge height of Model four is also higher than that of Model 2 after synkinematic sedimentation. Under the effect of synkinematic sedimentation, the development of these faults in the front of thrust wedges was promoted instead of new forward-vergent faults generated in the pre-wedge region ( Figure 6). The graphs of Figures 8A-C highlight the effect of synkinematic sedimentation in building a high angle and even supercritical wedge taper as well as in delaying the propagation of the deformation front into the foreland.

Analysis of Strain and Stress
In these four DEM models, the strain (volumetric strain and distortional strain) ( Figure 9) and stress (mean stress and maximum shear stress) ( Figure 10) characteristics of strata and faults in thrust wedge, pre-wedge, and undeformed regions are calculated after every 2 units of shortening.
The evolution of volumetric and distortional strain shows that: 1) At the beginning of new fault development, the volumetric and distortional strain are mainly concentrated in the shallow strata (Panels a3, b3, and d3 in Figures 9A,B,D). With the increase of model shortening, it gradually extends to the deep strata until it reaches the bottom detachment strata (Panels a4, b4, and d4 in Figures 9A,B,D). 2) In the thrust wedge, with fault F2 as the boundary, the volumetric strain between fault blocks in the hanging wall of fault F2 is mainly dominated by volumetric dilation (white to bright yellow), while between fault blocks in the footwall of fault F2, it is mainly dominated by volumetric contraction (light blue to white). 3) In the distortional strain field, along the shortening direction, the fault deformation is dominated by the top to the right sense of shear, showing the forward spreading fault, while the top to the left sense of shear is less, which causes the development of the backward spreading fault to be not obvious (Panels a3, b3, c3 and d3 in Figure 9). 4) synkinematic sedimentary loading has an obvious inhibitory effect on the evolution of new buried thrust faults in the prewedge region ( Figure 6 and Figure 7), causing almost no distortional strain to occur in this region ( Figures 9C,D c, d), relative to the other two models ( Figures 9A,B). However, sedimentary loading might have a promoting effect on these pre-existing faults in the front of the thrust wedge (such as the fault F3 in both Model three and Model 4, Figure 6 and Figure 7), with the increasing of these strain values along the fault plane, relative to other pre-existing faults F3 in Model 1 ( Figure 9A). It is also proved by the fault slip displacement of thrust fault F3 in Model three and Model 4 (Figure 7 c, d). 5) Compared with the distribution characteristics of distortional strain in these four models (Figure 9), the effect of erosion on the reactivation of preexisting faults (old faults F2 and F3 in the thrust wedge belts) could be presented by the widening of this distortional strain concentration zone along the fault plane. Moreover, the effect of erosion is more obvious for this old pre-existing fault near the moving plate wall (such as the fault F2 in Model two and Model 4) ( Figures 6B,D, 7B,D).
The mean stress and maximum shear stress reflect the overall characteristics of stress distribution and the instantaneous shear Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 659682 strength of the strata, respectively. In our models, stress analysis results show that: 1) The distribution of the high value of mean stress in each model is mainly concentrated at the bottom of these faults (Figure 9). Compared with Model 1, the high value of mean stress increased in different degrees in the other three models, especially in Model 3 with synkinematic sedimentation ( Figure 9C). This is consistent with the activity characteristics of thrust fault F3 in Model 3 ( Figure 7C).

Influence of Erosion and Sedimentary Loading on Activity of Thrust Faults in LMS
In active fold-thrust belts around the world, such as the Aconcagua fold-thrust belt in Argentina (Hilley et al., 2004), the central Andes fold-thrust belt in Bolivia (Horton, 1999;McQuarrie et al., 2008), the central Apennines fold-thrust belt in Italy (Scisciani and Montefalcone, 2006;Wu and McClay, 2011), and the Taiwan fold-thrust belt in China (Dahlen and Suppe 1988;Dadson et al., 2003;Wu and McClay, 2011), the roles of different synkinematic erosion and/or sedimentary loading, probably caused by long-term along-strike climatic variations (McQuarrie et al., 2008), on the generation and evolution of thrust faults in these fold-thrust belts (thrust wedges) have been interpreted and discussed by the Coulomb wedge model theory (e.g., Davis et al., 1983;Dahlen et al., 1984;Dahlen and Suppe 1988) and experimental simulations (e.g., Storti and Mcclay, 1995;Persson and Sokoutis, 2002;Simpson, 2006;Graveleau and Dominguez, 2008;Cruz et al., 2010;Cruz et al., 2011;Wu and McClay, 2011;Malavieille, 2011;Steer et al., 2014;Sun et al., 2016;Sun et al., 2021;Luo et al., 2021;Mao et al., 2021). Theoretical analysis and modeling results have all shown that syntectonic erosion reduced the number of major forwardvergent thrusts, and increased exhumation and thrust activities at the rear of the thrust wedge, resulted in out-of-sequence thrusting and fault reactivation in the wedge hinterland, and inhibited the forward propagation of the deformation front into the foreland. Meanwhile, the building of wedge taper by synkinematic sedimentation lowered thrust activities in the foreland relative to the hinterland. Similar results were obtained in our DEM experiments, such as promoting activity of thrust faults F2 and F3 in Model two and Model 4 ( Figures 6B,D, 7B,D, 9B,D), as well as inhibiting the forward propagation of thrust fault F4 in Model two and Model 4 ( Figures 6B,D, 7B,D, 9B,D), relative to Model one in this paper ( Figures 6A, 9A). These DEM results indicated that synkinematic erosion promoted the reactivation of the pre-existing faults in thrust wedges, and also have a negative influence on the formation and development of new incipient faults in the prewedges ( Figures 6B,D, 7B,D, 9B,D). Meanwhile, synkinematic sedimentary loading also has an effect on delaying the development of new incipient faults (such as thrust fault F4 in Model three and Model 4) in the pre-wedge region by promoting the development of pre-existing thrust fault(s) in the front of thrust wedge (such as the thrust fault F3 in Model three and Model 4) ( Figures 7C,D), causing the thrust wedge in a supercritical stage ( Figure 8A) with a narrow wedge length, relative to other models without synkinematic sedimentation ( Figure 8B).
The LMS fold-thrust belt is also an active fold-thrust belt with differential regional erosion in the interior (Gao et al., 2016;Tan et al., 2017) and sedimentation in the Sichuan Basin (Li Z. et al, 2018;Luo et al., 2021). The occurrence of the 2008 M W 7.9 Wenchuan earthquake has inspired a lot of thinking to the research of the LMS fold-thrust belt (Yin, 2010). One of these important questions is that the correlation between surface processes (such as erosion and sedimentary loading) and fault activity in these imbricate thrust systems of the LMS active foldthrust belt.
In this paper, according to modeling results of the activity of thrust faults F1, F2, F3, and F4 from DEM experiments, we built a comprehensive diagram to present the influences of erosion and sedimentary loading on fault activities in the LMS thrust wedge under the role of the compression stress ( Figure 11). We mainly discuss and analyze the characteristics of active faults in the front of the LMS fold-thrust belt, such as the WMF, YBF, PGF, and RFBT in the central LMS, and the JTF, WLF, SDF, and YAF (RFBT) in the southern LMS ( Figures 1B, 12). Here, we use the structural evolution of Model two and Model four to reflect the activity of faults in the southern and central LMS, respectively. Thrust fault F1 represents small faults that developed in shallow strata ( Figure 11) in the thrust wedge of the central and southern LMS. Thrust fault F2 and F3, which are large faults extending to the deep detachment layer (Figure 11), represent these preexisting thrust faults developing within and in the front of the thrust wedge, such as these interior faults WMF, JTF, WLF, YBF and frontal faults SDF, PGF in the central and southern LMS, respectively. Thrust fault F4 represents new incipient faults that developed as range front blind faults developed in the pre-wedge region, such as the buried faults RFBT and YAF in the front of the central and southern LMS.
According to the influence of erosion and sedimentation on the fault activity we described above, we can infer that the characteristics of erosion and sedimentation in the central and southern LMS have important influences on the activities of the WMF, YBF, PGF, and RFBT faults, and the WLF, SDF, and YAF faults, respectively. The activities of these faults are improved to some extent, which is consistent with the distribution characteristics of average erosion rates in the central and southern LMS (Figure 12). Perpendicular to the strike of LMS, different average erosion rates (AER C1 , AER C2 , and AER C3 in the central LMS and AER S1 and AER S2 in the southern LMS, Figure 12) promoted the reactivation of these faults to different extents under the effect of regional compression stress. In our models, the activity of thrust fault F1 was ended early and kept inactive in the process of the development of thrust wedge, whether it was under the effect of synkinematic erosion or not. It infers that these faults developed in the shallow strata of fault blocks might be inactive in the process of the development of the LMS fold-thrust belt. However, different structural evolution occurred on thrust faults F2 and F3, which has a large amount of slip displacement during the evolution of thrust wedge. It indicates that the YBF and PGF faults in the central LMS, and the WLF and SDF faults in the southern LMS are major active faults, especially under the effect of synkinematic erosion (Figures 6, 9, 10). It was also confirmed by Tan et al. (2017) based on apatite fission track (AFT) ages in the Xuelongbao Massif and the Pengguan Massif, which indicate Quaternary thrust activity of the WMF. Moreover, according to the activity of the thrust fault F3 (Figure 11), we can infer that PGF in the central LMS and SDF in the southern LMS has the characteristics of a large amount of slip displacement with persistent activity at present, especially under the role of synkinematic sedimentation deposited in the range front ( Figures 7C,D). Gao et al. (2016) indicated that the YBF in the Central LMS and the SDF in the Southern LMS act as major tectonic boundaries separating areas experiencing rapid uplift from slow uplift based on the results of the geomorphic analysis. It is consistent with experimental results in this paper.   Gao et al., 2016). AER data (AER C1 , AER C2 , and AER C3 are decreasing successively in the central LMS; AER S1 is bigger than AER S2 in the southern LMS) from Figure 2 and Figure 3 in this paper. WMF, Wenchuan-Miaoxian Fault; YBF, Yingxiu-Beichuan Fault, PGF, Pengxian-Guanxian Fault; JTF, Jiatang Fault; WLF, Wulong Fault; SDF, Shuangshi-Dachuan Fault; YAF, Ya'an Fault. Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 659682 In addition, the experimental results show that the sedimentary loading in the range front (such as the Chengdu plain) also has an important influence on the structural evolution and fault activity of the LMS fold-thrust belt ( Figure 11). In the range front of the southern and central LMS fold-thrust belt, synkinematic sedimentary strata were deposited with different thickness (Figures 1C, 2B, 3B). Our experimental results show that the synkinematic sedimentary strata can prevent the development of new leading faults thrusting forward in the pre-wedge region, such as the thrust fault F4 ( Figures 6B,C, 7B,C). Hence, the differential sedimentary distribution may result in the differential evolution of thrust wedge in the central and southern LMS. We will focus on this part in Differential Evolution of Thrust Wedge in the Southern and Central LMS. At the same time, synkinematic sedimentation also has a good blocking effect on the stress propagation and strain convergence on the fault planes, so that the highest stress on the fault planes basically converge at the root of the main active faults inside the thrust wedge ( Figures 10C, 11), forming a good seismogenic environment in this area. Hence, relative to northern LMS, the central LMS and southern LMS is more conducive to the occurrence of earthquakes, because of the growth of the Chengdu plain, a structural basin filled with synkinematic sedimentation. Actually, the central LMS and southern LMS is a concentrated area of epicenter distribution (Figures 2, 3) (Liu-Zeng et al., 2009;Xu et al., 2009Xu et al., , 2013Li et al., 2010;Wang et al., 2014).

Differential Evolution of Thrust Wedge in the Southern and Central LMS
As we discussed in Influence of Erosion and Sedimentary Loading on Activity of Thrust Faults in LMS, the synkinematic sedimentary loading can prevent the development of new leading faults thrusting forward in the pre-wedge region by promoting the development of pre-existing faults in the front of thrust wedge ( Figures 6B,D, 7). If a certain thickness of synkinematic sedimentation was deposited in front of the thrust wedge, the width of the wedge decreases, while the height of the wedge increases after the same shortening displacement ( Figures 8B,C). For example, under the same erosion rate, different syntectonic sedimentary thickness has an important influence on the overall development of the thrust wedge, which causes the occurrence of the geomorphic characteristics of recess and salient in the front of thrust wedge during the compression process ( Figure 13).
These differences are in good correspondence with the current tectonic and geomorphic characteristics of the central and southern LMS (Figures 12, 14). From the central LMS to the southern LMS, the Cenozoic depocenter is located in the central LMS, which is at least 200 m thicker than that in the southern LMS ( Figure 1C). This kind of differential deposition may lead to different characteristics of fault activity and uplift in the interior thrust wedge and prewedge region in the central and southern LMS, such as the differential fault activity of the PGF and SDF faults, and the RFBT fault. While at the same time, the differential evolution of these faults was likely to result in differential geomorphic evolution (Tan et al., 2018), such as the occurrence of Dujiangyun recess (Liu Y. et al., 2020). Liu Y. et al. (2020) has already highlighted these geomorphic characteristics, such as the formation of salients and recesses in the front of the LMS fold-thrust belt. Moreover, the results of analog experiments from Luo et al. (2021) have also shown that the deposition and distribution of the wedge-shaped quartz sand (e.g., alluvial fan in nature), caused by regional erosion (e.g., fluvial incision in nature), played an important role in controlling along-strike variations in thrust-belt deformation, causing local variations of thrust-belt geometry. Therefore, except for the erosion, the differential sedimentary loading may also be one of the important factors for the differential evolution along the strike of LMS fold-thrust belt, such as the differential FIGURE 13 | The DEM results in this paper (A) and the conceptual model modified from Liu et al., 2020 (B) present the occurrence of recess and salient in the front of thrust wedges.

Limitations of DEM Experiments
Our DEM experimental set-ups include some simplifications that ensure the feasibility of the study and allow the investigated parameters to be highlighted. But we need to keep the limitations in mind when applying the results to natural examples. 1) The structural evolution of the LMS fold-thrust belt is more complex than a single wedge, e.g., Sun et al. (2018) highlighted the role of pre-existing structures in shaping the Cenozoic tectonics of the LMS and Sun et al. (2016) suggested that across-strike preexisting topographic relief controls the deformation partitioning within the LMS fold-thrust belt. In this paper, pre-existing structures in the western LMS fold-thrust belt were not implemented in the current discrete element models. In addition, only several main faults in the LMS fold-thrust belt were simulated and analyzed in extremely simplified models and parameters of particles. 2) As found out by Hubbard et al. (2010), two taper wedges (topographic slopes) could be defined from the hinterland to foreland of the LMS fold-thrust belt: A steep surface slope caused by a slightly strong basal detachment in the front of the LMS fold-thrust belt (e.g., Chen and Wilson, 1996;Li et al., 2014), and a shallow surface slope caused by a shallow or middle level, weak detachment (e.g. the Triassic and Cambrian salts, and the Silurian shale) in the Sichuan basin (e.g., Cai and Liu, 1997;Jia et al., 2006;Tang et al., 2008;Li Z.-G. et al., 2013;Sun et al., 2021). Moreover, Sun et al. (2021) have presented the role of the middle level detachment in shaping the current western LMS fold-thrust belt by analogue models. In our models, we only focused on the formation of the thrust wedge (the front of the LMS fold-thrust belt) by considering the deep slightly strong basal detachment (contrasting with the salt detachment). The structural evolution in the Sichuan basin with a shallow or middle level, relatively weak detachment (evaporite or -shale) was not modeled during the shortening. 3) To better understand the role of the single parameter in the development of the thrust wedge, the erosion and the sedimentary loading processes were modeled separately in model two and model 3. Although in model 4, the erosion and the sedimentary loading processes were modeled together, thees were still not very natural. For example, they are not continuous in the models, which is more or less different from nature. Also, erosion was achieved by removing everything above a certain height, and no erosion occurred on the wedge below this cap. Sediments in the foredeep of a foreland basin usually had a wedge shape in nature, but in our models a constant thickness was used across the basin.

CONCLUSION
Four groups of DEM models were set up to simulate and analyze the influence of regional erosion and sedimentary loading on the formation and spatial-temporal evolution of faults in the southern and central LMS active fold-thrust belt. The interior characteristics of faults in the southern and central LMS foldthrust belt were evaluated during the interaction of tectonic processes and surface processes according to the stress-strain analysis from DEM results. Model results demonstrated that: (1) Synkinematic erosion promoted the reactivation of the preexisting faults in thrust wedges and had a negative influence Frontiers in Earth Science | www.frontiersin.org on the formation and development of new incipient faults in the pre-wedges. Synkinematic sedimentary loading also delayed the development of new incipient faults in the prewedge regions by promoting the development of thrust faults in the front of thrust wedges, causing these thrust wedges in supercritical stages with narrow wedge lengths, relative to other wedges without synkinematic sedimentation.
(2) The characteristics of erosion and sedimentation have important influences on the activities of large faults which are extended into the deep detachment layer, such as the WMF, YBF, and PGF faults, and the WLF, SDF faults in the central and southern LMS, respectively. Perpendicular to the strike of LMS, different average erosion rates promoted the uplift of these faults to different extents. Moreover, due to the sedimentary strata in the range front or foreland basin, the sedimentary loading plays a role in preventing or delaying the reactivation of the underlying pre-existing faults or the generation of new incipient faults, inferring that the evolution of the RFBT or YAF in the central or southern LMS is inhibited and delayed because of the growth of the Chengdu plain at present. (3) Besides differential erosion, the differential sedimentary loading may also be one of the important factors for the differential evolution along the strike of the LMS fold-thrust belt. This kind of differential deposition may lead to different characteristics of fault activity and uplift in the interior thrust wedge and pre-wedge region in the central and southern LMS. (4) Synkinematic sedimentation has a good blocking effect on the stress propagation and strain convergence on the fault planes, so that the highest stress on the fault planes converges at the root of the main active fault inside the thrust wedge, forming a good seismogenic environment in this area. Therefore, relative to the northern LMS, the central LMS and southern LMS is more conducive to the occurrence of earthquakes.

ACKNOWLEDGMENTS
We would like to thank Julia Morgan and Thomas Fournier for generously sharing their post-processing scripts and algorithms, which have been used to process and display the model outputs presented here. We would like to express our gratitude to these reviewers for their careful and helpful comments, which significantly improved the manuscript.