Stepwise Widening of the Central Andes—The Role of the Lower Crust

The outward growth of many orogens occurs through pulsed migration of the strain front. During Cenozoic shortening of the central Andes, the strain front abruptly migrated ∼400 km eastward across the Altiplano, isolating the plateau interior from major deformation. In contrast to the traditional critical wedge model that focuses on shallow crust deformation, our lithosphere‐scale numerical models show the ductile lower crust plays a key role in the dynamics of the orogen. As the orogen grows upward, the lower crust flows outward and causes stepwise expansion of the orogen through landward‐migrating deformation. The step length of the strain‐front migration increases with decreasing lower crustal strength. Topographically, a weaker crust promotes the formation of a smooth plateau and underthrusting of stronger foreland lithosphere. Our results indicate that the seemly unrelated physiography, deep‐lithosphere structure, and plateau‐wide strain migration in the Altiplano may all be the result of a weak lower crust.

crust may be ductile, which allows it flow laterally under stresses created by topographic loads, strongly influencing the orogen and wedge structure (Beaumont et al., 2006;Braun & Yamato, 2010;Clark & Royden, 2000;Vanderhaeghe et al., 2003). It is unclear how this ductile crustal layer may influence the kinematics of orogenic growth, mountain belt physiology and deep lithosphere structure.
Here, we focus on the central Andes, which is the largest plateau on Earth formed through ocean-continent convergence (Figure 1a). The Altiplano of southern Bolivian (ca. 18°-23°S) is the widest region (∼700 km) in the central Andes. It formed through ca. 300-400 km of crustal shortening, mainly during last 40 Ma (Eichelberger et al., 2015). At present, the Altiplano has a surface elevation of ca. 4 km, with low relief in the plateau interior. The topography gradually decreases eastward from the high plateau to the low Chaco plain over a distance of ∼300 km.
Reconstruction of the foreland basin system in the Altiplano indicates that the orogen was narrow, with shortening limited to the Precordillera between the Late Cretaceous and mid-Paleocene (Figure 1b;DeCelles & Horton, 2003). After ∼25 Myr of shortening in the Precordillera, the Eocene-to-Oligocene stratigraphy indicates that the forebulge deposits abruptly changed to foredeep deposits in the Eastern Cordillera, which suggests a rapid eastward passage of the foreland basin (DeCelles & Horton, 2003;Horton, 1998). Thermochronological data from the Eastern Cordillera also show rapid cooling and exhumation, which was possibly caused by thrust-driven rock uplift in the Eocene (Ege et al., 2007). These independent observations suggest that the retroarc thrust belt abruptly migrated from the Precordillera to the Eastern Cordillera in the middle-to-late Eocene (e.g., Elger et al., 2005;Horton, 2005;McQuarrie, 2002). This ∼400 km eastward jump of the retroarc strain loci isolated the Altiplano interior from deformation during the orogenic widening event (McQuarrie et al., 2005).  (Beck & Zandt, 2002;Dorbath & Granet, 1996;Myer et al., 1998;Metcalf & Kapp, 2015 and references therein). PrC = Precordillera; WC = Western Cordillera; EC = Eastern Cordillera; IA = Interandean zone; SA = Subandes; SB = Santa Barbara Ranges; SP = Sierras Pampeanas. (b) Migration of retroarc orogenic strain loci at 21°S (Quiero et al., 2022;Stalder et al., 2020). Orange boxes represent compression, dotted black boxes indicate extension. Blue dashed arrows show the migration of retroarc strain front. (c) Schematic cross section for the Altiplano (Beck & Zandt, 2002;Whitman et al., 1996). The zoom in plot shows the P wave velocities from Myers et al. (1998), where the high velocities beneath the eastern Altiplano plateau are interpreted to be the underthurst Brazilian craton.
After this major orogenic front migration, the fold-thrust deformation stalled in the Eastern Cordillera from ca. 40 to 36 Ma and developed a steeply dipping, bivergent thrust (McQuarrie & DeCelles, 2001;Müller et al., 2002). Starting at ca. 35 Ma, deformation in the Eastern Cordillera gradually ceased, and foredeep deposition started to advance eastward (DeCelles & Horton, 2003;Horton, 2005). The thrust belt propagated steadily across the Interandean zone from ca. 25 to 10 Ma and the Subandean zone from ca. 10 Ma to present (Horton & DeCelles, 1997;Uba et al., 2006). These kinematics created a broad thin-skinned eastward-vergent thrust belt in the foreland (Figure 1c; Kley, 1996). The general retro-wedge kinematics since the Late Cretaceous are roughly consistent along strike between 15°S and 23°S .
The stepwise widening of the Altiplano through major jumps in strain loci has been attributed to different factors. One proposal is that this is controlled by inherited structures from ancient rifts in the Eastern Cordillera, which are weak enough to be the locus of the new orogenic/strain front (Perez et al., 2016;Sempere et al., 2002;Viramonte et al., 1999). Other studies interpret the large migration of the Altiplano strain front to be caused by dynamic events, such as a low-angle subducting slab (Martinod et al., 2020) or episodic lithosphere removal (DeCelles et al., 2009).
Here, we investigate an alternate factor that may control orogenic growth-a ductile lower crust. Continental lithosphere can have a stratified rheological structure, with a weak lower crust sandwiched between a strong upper-mid crust and strong mantle lithosphere (Kohlstedt et al., 1995;Ranalli & Murphy, 1987). Studies show that a weak lower crust can influence intracontinental deformation (e.g., Bendick & Flesch, 2007;Wang & Currie, 2017). Previous work has investigated the conditions for crustal flow (Gerbault & Willingshofer, 2004;McKenzie et al., 2000), the redistribution of crustal materials by viscous flow to account for discrepancies between shortening rate and crustal thickness (Gerbault et al., 2005;Husson & Sempere, 2003;Ouimet & Cook, 2010;Yang et al., 2003), and the influence of crustal flow on surface morphology (Clark et al., 2005;Clark & Royden, 2000;Cook & Royden, 2008). However, the role of the ductile lower crust in the nucleation of brittle shear zones remains unclear. In this study, we explore how the lower crustal strength controls the first-order kinematics of the orogenic strain front and its effects on surface geomorphology and lithosphere deformation.

Numerical Models
We designed two sets of thermal-mechanical coupled numerical models using the finite element code SOPALE (Beaumont et al., 2006) (details in Supporting Information S1). The 2D models examine the shortening of continental lithosphere using materials that have frictional-plastic and viscous rheologies. Plate shortening is introduced by imposing an uniform influx velocity at the right boundary of lithosphere at 8 mm/yr (Figures 2a and 3a), comparable to the Cenozoic shortening rate in the central Andes . We first present simplified models with a constant viscosity in the lower crust to investigate the relationship between lower crustal strength and development of the orogenic strain front. Then more realistic models are presented that are based on the central Andean structure and include a non-Newtonian lower crust and a thick craton in the foreland. These enable us to examine the effects of strain-dependent lower crust deformation on orogen evolution and the dynamics at the boundary between the orogen and craton.

Simplified Models
In the simplified models, the lithosphere is divided into a 200-km-wide strong forearc block and a 1000-km-wide weak block (Figure 2a). The plastic cohesion and viscosity of the forearc block are 10 times higher than the weak block, and thus deformation is confined to the weak region, consistent with the observations that shortening mainly occurred to the east of the Andean forearc. The 100-km-thick lithosphere includes 24 km upper crust and 16 km lower crust. A 100-km-wide region with elevated temperatures is placed near the west margin of the weak block to simulate the volcanic arc. A constant viscosity is used for lower crust in the weak block; values ranging from 10 19 to 10 21 Pa s are tested. The upper crust, mantle lithosphere and sublithospheric mantle have non-Newtonian rheologies based on laboratory studies (Table S1 in Supporting Information S1). Field and laboratory observations suggest that deformation-induced grain size reduction, kinematic metamorphic reactions, and/or fluid infiltration can cause strain softening and weakening in the brittle and ductile regimes (e.g., Behn et al., 2002;Bos & Spiers, 2002;Jin et al., 1998;Karato et al., 1986). Strain softening and weakening are included in the models through linearly decreasing the viscosity of the upper crust and mantle lithosphere by a factor of 10 A 100-wide region with elevated temperature is placed at the western side of weak block to represent the volcanic arc area (dashed red box). Model runs have two stages. Models first undergo isostatic adjustment due to subsurface density variations. Then plate shortening is imposed by pushing the right-side boundary of lithosphere at 8 mm/yr. The left-side boundary of lithosphere is fixed. WQ = wet quartzite (Gleason & Tullis, 1995); MG = mafic granulite (Wang et al., 2012); ML = mantle lithosphere; DO = dry olivine (Hirth & Kohlstedt, 2003); WO = wet olivine (Karato & Wu, 1993); temp. = temperature. (b) Surface topography and cumulative shear strain in the lithosphere for models where the lower crust (LC) has constant viscosities of 10 19 , 10 20 and 10 21 Pa s, respectively. All times are since the start of plate shortening. The zoom in plots show the horizontal deviatoric stress (upper panel; positive values denote compression and negative values denote extension) and the crustal horizontal velocity field (lower panel) shortly before the development of the nascent shear zone. After the new strain front forms at the orogenic tip, the compressional stress decreases. Profiles of horizontal velocity are also shown at the locations indicated by the vertical lines. (c) The model surface topographies at the time that the new strain fronts develop. Different color lines represent models with different viscosities in the lower crust, ranging from 10 19 to 10 21 Pa s. Numbers in brackets are the times when the new shear zones nucleate. and reducing the effective angle of friction ϕ eff from 8° to 2° as the accumulated strain ε increases from 0.1 to 1 (Huismans & Beaumont, 2003;Warren et al., 2008).
Crustal shortening is initially localized in the thermally weakened arc region, resulting in crustal thickening and surface uplift. The high topography creates high pressures in the underlying crust, causing the lower crust to flow eastward ( Figure 2b). As a result of the outward crustal flow, the slope of the orogenic flank shallows and the orogenic wedge grows landward. The surface topography, resulting from both plate shortening and lower crustal flow, redistributes the stresses in crust. After the arc surface uplifts to 1.8-2 km, the horizontal stress changes in the arc region from compression to extension (Figure 2b). Meanwhile, the toe of orogenic wedge experiences compression. When the stress exceeds the crustal strength, the rock fails and a nascent strain zone develops at the wedge toe. The location of the new strain front shows a strong correlation with the lower crustal strength. For a lower crustal viscosity of 10 19 Pa s, the new strain front nucleates ∼550 km away from the original strain loci at the arc at ∼11.6 Myr (Figure 2b). For a lower crust viscosity of 10 20 or 10 21 Pa s (Figure 2b), the eastward flow in the lower crust is limited by the increased viscosity. Consequently, the orogenic flank steepens and the distance between the arc and new strain front reduces to ∼200 and ∼150 km, respectively. The nucleating times of the new shear zones are also reduced to ∼8.8 and ∼6.8 Myr respectively, due to the limited outward crustal flux. Additional lower crustal viscosities (5 × 10 19 and 5 × 10 20 Pa s) are tested and the distance and nucleation time of the new strain front systematically become larger as the lower crust becomes weaker (Figure 2c).
Our models demonstrate additional lithosphere dynamics. Below the high-temperature region, shortening induces gravitational foundering of the lowermost mantle lithosphere, which is negatively buoyant as it is cooler than the underlying mantle (Figure 2b). The small-scale drips entrain the adjacent lithosphere and induce asthenospheric upwelling, but there is little effect on the surface elevation. This is because the downwelling lithosphere has a small negative thermal buoyancy and therefore the viscous stresses on the overlying material are low (Conrad & Molnar, 1999). Furthermore, the weak lower crust inhibits the upward transmission of mantle stresses (Wang & Currie, 2017).

Central-Andean-Type Models
In this series of models, all materials have non-Newtonian and temperature-dependent viscosities. The lower crust uses the laboratory-derived rheology of mafic granulite (Wang et al., 2012); all other material parameters are identical to those in the simplified models. The weak block is initially 500 km wide and a 200-km-thick craton lithosphere is placed at the eastern side of the weak block (Figure 3a), representing the thick Brazilian foreland lithosphere. The Moho temperature is ∼740°C in the weak block and ∼410°C in the craton. As the rheology and density parameters used for craton are identical as those in the weak block, the relatively low temperatures make the craton stronger and denser.
In Model 1, the lower crust is 5 times weaker than the laboratory mafic granulite (Table S1 in Supporting Information S1), simulating the strength of a more felsic composition or a higher degree of hydration than the laboratory sample, or a hotter thermal structure. In this model, initial crustal shortening occurs at the arc. Beneath the uplifted arc, the upper crust moves westward at ∼5 mm/yr due to plate shortening. However, the lower crust flows in the opposite direction at a rate of up to ∼10 mm/yr at ∼5 Myr (Figure 3b). To the east of arc, the lower-crustal flow velocity gradually decreases, and the upper crust at the eastern margin of weak block experiences compression. After ∼120 km of crustal shortening, an east-verging thrust belt develops at the eastern margin of the weak block. This new shear zone is ∼350 km away from the original strain loci at the arc.
The crustal strength evolves as the strain rate varies during plate shortening. Beyond the hot arc, the lower crustal viscosity is >5 × 10 20 Pa s in the weak block at the beginning of shortening and decreases over time ( Figure  S1 in Supporting Information S1). As the new strain front develops at the eastern flank of weak block at ∼15 Myr, the lowermost crust across the weak block is ∼10 19 Pa s (Figure 3b and Figure S1 in Supporting Information S1). The weak lower crust decouples the upper crust and mantle lithosphere. Compared to the overlying lower crust, the mantle lithosphere beneath the weak block continues moving westward. Small-scale drips gradually develop at the base of lithosphere at ∼6 Myr. The drips and the base of residual lithosphere are entrained eastward by the edge-driven flow created by the craton (Figure S2 in Supporting Information S1 and Movie S1). As the weak-block crust thickens, the underlying mantle lithosphere does not significantly thicken, owing to the drips. A comparison of topography before and after lithosphere removal shows that the small-scale drips do not result in significant surface uplift, due to the small "available buoyancy" of the downwelling lithosphere (Conrad & Molnar, 1999) and because the weak lower crust inhibits the upward transfer of mantle stresses (Wang & Currie, 2017). Therefore, variations in surface topography are mainly controlled by variations in crustal thickness and the kinematics of the crustal shear zones.
After the development of the east-verging thrust fault, the orogenic deformation stalls at the transition between the orogen and craton. Instead of continuing to propagate cratonward, the deformation moves back into the orogen and develops a west-verging shear fault at ∼18 Myr (Movie S1). The plateau uplifts until the gravitational weight becomes sufficient to cause the lower crust to flow into the craton. A thin area with viscosity 10 20 -10 21 Pa s gradually develops in the lowermost crust of the craton. At shallow depths, a series of east-verging thrust faults migrates into the craton crust starting at ∼24 Myr. The weakened lower crust allows the craton mantle lithosphere to decouple from the overlying crust, and it begins to underthrust westward at ∼20 Myr. At 40 Myr (320 km plate shortening, comparable to the estimated total shortening in the Altiplano; Eichelberger et al., 2015), a ∼100-km-wide craton mantle lithosphere underlies the crust of weak block.

Variations in Lower-Crustal Strength and Weak Block Width
Model 1 shows that orogen deformation is influenced by ductile deformation of the lower crust in the weak block and the presence of a strong craton that can obstruct the crustal flow. We have tested a range of lower-crustal strengths and widths for the weak block. Models show two main behaviors (Figure 4a): (a) When the lower crust is weak and the weak block is narrow, crustal flow expands across the weak block but is inhibited by the craton, and the new shear zone nucleates at the boundary of the craton (e.g., Model 1). (b) When the crust is strong or the weak block is wide, the lateral expansion of crustal flow only occurs in a portion of the weak block before the new strain zone develops in the interior of the weak block. For example, when the lower crust is 10 times stronger than that in the Model 1, the new strain front initiates ∼120 km away from the arc at ∼16 Myr (Model 2; Figure 4b and Movie S2). After this, each migration of the orogenic strain front is less than 50 km. After 240 km of shortening in Model 2 (comparable to the shortening amount in the Puna region of the central Andes; Anderson et al., 2017 and references therein), the interior of the plateau is deformed by several bivergent shear zones (Figure 4b).
Compared to the widespread weak lower crust in Model 1, the weak crust (≤10 20 Pa s) in Model 2 is spatially limited to the lowermost crust and develops slowly cratonward (Figure 4b and Figure S3 in Supporting Information S1). The new shear zone nucleates within the weak block. At 30 Myr, the lower crust viscosity beneath the craton remains >10 22 Pa s, and consequently, the crust and mantle lithosphere are strongly coupled, resulting in the lack of underthrusting of foreland lithosphere (Figure 4b).

Variations in Shortening Rate, Brittle Strength and Strain Softening
Model sensitivities to the plate shortening rate, brittle strength, and strain softening parameters are tested in additional models where the lower crust is either viscously weak (Model 1) or strong (Model 2). With a weak lower crust, reducing the shortening rate provides more time for crustal flow, which widens the mountain range and smooths the ground surface (Model 1-a; Figure S4 in Supporting Information S1). If the brittle crust is stronger (i.e., higher initial friction angle), rock failure is inhibited and deformation tends to localize in the primary strain region (Model  1-b). In both cases, if the lower crust is stronger, the migration lengths of strain loci are reduced (Models 2-a and 2-b; Figure S5 in Supporting Information S1). If there is no frictional-plastic softening (Model 1-c) or viscous weakening (Model 1-d), the lack of positive feedback between strain accumulation and lithospheric weakening inhibits localized deformation and results in numerous shear zones ( Figure S4 in Supporting Information S1). This is consistent with previous work that shows that strain softening promotes the development of focused shear zones (e.g., Huismans & Beaumont, 2003). When there is no strain softening in mantle lithosphere (Model 1-e), the dynamics of strain loci are similar to those in Model 1 ( Figure S4 in Supporting Information S1). In this case, a stronger lower crust also decreases the migration distance of strain loci (Model 2-e; Figure S5 in Supporting Information S1), consistent with the relationship between lower crustal strength and kinematics of strain front shown in our other models.

Discussion and Conclusions
The episodic, stepwise migration of the retroarc orogenic front is an unsolved question in continental tectonics. In the Altiplano region of the central Andes, stratigraphic, thermochronologic and structural data suggest that shortening was accompanied by a plateau-wide eastward jump of the orogenic front in the Eocene (McQuarrie et al., 2005). In this study, we explore the role of a laminated crustal rheology in orogenic growth. Our models show that the strength of lower crust not only influences the surface morphology (Clark & Royden, 2000), but also controls the location of the strain front. As plate shortening creates a localized area of thickened crust, the gravitational weight of the high topography causes the ductile lower crust to flow outward, and a new strain zone  Wang et al. (2012). λ is the ratio of the migration distance of the strain loci to the width of shortened weak block when the new strain loci nucleates. Circles denote models in which the strain loci migrate across the weak block to the boundary between weak block and craton. Triangles denote models in which the strain loci forms within the interior of the weak block. (b) Model 2 results. The lower crust is 10 times stronger than that in Model 1. Model structure and boundary conditions are the same as those in Model 1. Left column shows the model evolution and inset plots shows the horizontal velocity. Vectors show the particle motion. Right column shows the cumulative shear strain in the lithosphere and inset plots show the horizontal deviatoric stress and profiles of lithospheric viscosities. nucleates due to the compressional stresses at the toe of orogenic wedge. The distance of the new strain front increases with decreasing lower crustal strength. If a cold craton blocks the crustal flow, the compressive stress localizes at the boundary between the orogen and craton.
Our models indicate that the formation of east-verging thrust faults in the Eastern Cordillera during the Eocene (Ege et al., 2007) may be related to eastward lower crust flow beneath the Altiplano that was blocked by the strong Brazilian craton. The models show that when the surface load of the plateau is insufficient to drive lower crustal flow into the craton, crustal deformation moves backward into the plateau, developing a backthrust shear zone at the eastern side of the weak block (Figure 3b at 20 Myr). This may explain the out-of-phase backthrust zone that developed in the western Eastern Cordillera following the east-verging thrust belt to the east (McQuarrie & DeCelles, 2001). Thus, we propose the bivergent thrust structure reflects the stalling of plate shortening and lower crustal flow at the boundary between the plateau and craton.
In the models, the continued plateau shortening and uplift lead to the development of deformation that migrates into the cratonic crust. At the time that crustal deformation starts to migrate into the craton region, the eastern margin of the weak block has uplifted to ∼4 km (model time at 30 Myr in Figure 3b). This is consistent with a paleoelevation study that suggests the Eastern Cordillera of the Altiplano was close to its present-day elevation during the late Miocene when the fold-thrust deformation migrated into the Subandes (Leier et al., 2013). In the models, the base of cratonic crust weakens over time (<10 21 Pa s) and the weakening area gradually expands landward. This produces a series of east-verging shear zones within the craton that are closely spaced. Therefore, the observed kinematic change from the plateau-wide jump in the Altiplano to the progressive migration of deformation in the Interandea/Subandean zones (McQuarrie et al., 2005) may be caused by the differences in crustal temperature and strength between plateau and craton. In the deep lithosphere, the model crustal thickness is 60-75 km beneath the weak block at 40 Myr (320 km of shortening), with the thickest crust located at the eastern margin. To the east, the crust gradually thins to 40 km in the craton region. The model crustal thicknesses are comparable to the present-day seismically determined crustal thickness for the Altiplano . Furthermore, the models show that the weak lower crust in the craton functions as a crust-mantle décollement that allows the cratonic mantle lithosphere to decouple from the crust and underthrust the plateau. This is consistent with seismic results that indicate the Brazilian craton has underthrust westward to ∼65.5°W in southern Bolivia (Figure 1c; Beck & Zandt, 2002 and references therein).
Our models may also explain the difference in deformation between the Altiplano and Puna regions of the central Andes. The Puna is characterized by a rugged and fragmented hinterland, where deformation progressed in smaller steps. In our models, when the hinterland lower crust is relatively strong, ductile flow is confined to a narrower region and orogenic jumps are more closely spaced. This creates a rough plateau surface with steep orogenic flanks (e.g., Model 2), as seen in the Puna. In contrast, the weak lower crust in Model 1 allows for a relatively smooth surface in the plateau interior and a gentle flank where the surface elevation gradually decreases from 4 to 5 km to the low craton over a distance of ∼350 km (Figure 3b), comparable to the Altiplano topography. A stronger Puna crust may also inhibit underthrusting of the foreland lithosphere, which is consistent with seismic observations beneath the Puna (Gao et al., 2021).
This study focuses on the role of lower crust strength and provides a first-order description of the lithosphere deformation during orogenic shortening in two-dimensional models. We note that in three-dimensions the crust can flow along the strike of the orogen (e.g., Gerbault et al., 2005), which would reduce the velocity of across-strike crustal flow. To achieve the crustal flux and surface relief observed in our models, a higher shortening rate or more shortening time may be required to nucleate new shear zones. We also acknowledge that Andean evolution involves multiple processes that are not included in our models. Lateral heterogeneity in lithospheric strength (e.g., ancient rift zones; Sempere et al., 2002;Perez et al., 2016), surface process (e.g., erosion; Horton, 1999), subduction (Faccenna et al., 2017;Martinod et al., 2020), crustal melting and weakening (de Silva & Kay, 2018), changes in plate convergence (Quiero et al., 2022;Somoza, 1998) and removal of eclogitized lithosphere (Wang et al., 2021) also contribute to the detailed Andean evolution. Additionally, the thick sediments in the Subandean zone, which are not included in our models, may help to develop a thin-skinned thrust-belt system with a shallow detachment in the Altiplano Babeyko et al., 2006;Ibarra et al., 2019). More detailed models that assess the effects of these factors on the orogen deformation should be investigated in the future.
Despite the simplifications, our models show that ductile flow of lower crust significantly modifies the dynamics of orogen shortening and may explain the stepwise widening of orogens. A weak lower crust facilitates the migration of retroarc thrusting into the continental interior, and lower crustal strength controls the spacing and timing of the landward migration of deformation. This mechanism can explain the punctuated migration of deformation in the Altiplano region, and based on shared physical mechanisms, the results also provide insights into how other orogens grow and widen in nature.

Data Availability Statement
The data used in this paper are from numerical models with the model parameters shown in Tables S1 and S2 in Supporting Information S1. Model parameters and output data can be found in Zenodo repository (https://doi. org/10.5281/zenodo.7804452). Details of the model software are available at http://geodynamics.oceanography. dal.ca/sopale_nested.html. Figures were made with GMT v6.1.0 (https://www.generic-mapping-tools.org).