Hindered Trench Migration Due To Slab Steepening Controls the Formation of the Central Andes

The formation of the Central Andes dates back to ∼50 Ma, but its most pronounced episode, including the growth of the Altiplano‐Puna Plateau and pulsatile tectonic shortening phases, occurred within the last 25 Ma. The reason for this evolution remains unexplained. Using geodynamic numerical modeling we infer that the primary cause of the pulses of tectonic shortening and growth of the Central Andes is the changing geometry of the subducted Nazca plate, and particularly the steepening of the mid‐mantle slab segment which results in a slowing down of the trench retreat and subsequent increase in shortening of the advancing South America plate. This steepening first happens after the end of the flat slab episode at ∼25 Ma, and later during the buckling and stagnation of the slab in the mantle transition zone. Processes that mechanically weaken the lithosphere of the South America plate, as suggested in previous studies, enhance the intensity of the shortening events. These processes include delamination of the mantle lithosphere and weakening of foreland sediments. Our new modeling results are consistent with the timing and amplitude of the deformation from geological data in the Central Andes at the Altiplano latitude.

Despite the multitude of proposed shortening mechanisms, none adequately explains the evolution and variability of the deformation in the Central Andes during the last ∼35 Ma ( Figure 1c). However, the resolution of the shortening rate compilation (<5 Ma) from Oncken et al. (2006) offers a solid base to investigate this problem through geodynamic models. Although the data may carry intrinsic uncertainties from using different measurement methods, it shows a systematic consistency in shortening amplitudes across time and latitude.
Shortening rates along the Altiplano section at 21°S are the most temporally resolved and suggest four different phases of deformation in the last ∼50 Ma (Figure 1c). Between ∼50 and 33 Ma (Phase 1), the shortening rate linearly increased from 0 to ∼3.5 mm/yr before escalating at ∼33 Ma to ∼8 mm/yr. From ∼33 to 15 Ma (Phase 2), the shortening rate stagnates in a range between ∼4 and ∼7 mm/yr. From ∼15 to 7 Ma (Phase 3) the shortening rate pulsed, reaching a maximum of ∼11 mm/yr before dropping to ∼5 mm/yr. Following this, a second pulse occurred from ∼7 Ma to present (Phase 4) that reached a maximum value of ∼16 mm/yr before dropping to the ∼8 mm/yr seen from present day GPS velocities (Bevis et al., 2001;Klotz et al., 2006).
Utilizing high-resolution geodynamic models, with buoyancy-driven subduction, and validating them through geological shortening data from the Central Andes, this study proposes a new mechanism associated with slab steepening due to buckling events. This mechanism provides an explanation for the variability of the shortening rate. The models additionally address the decline in deformation intensity between 7 and 4 Ma to present-day levels. Our results suggest that a complex interaction between the oceanic and continental plates controls the timing and variability of the deformation in the Central Andes since the Oligocene.

Governing Equations
We use the geodynamic finite element code ASPECT (Advanced Solver for Problems in Earth's ConvecTion, Bangerth et al., 2021;Heister et al., 2017;Kronbichler et al., 2012;Rose et al., 2017) to setup a 2D subduction model. The model solves three conservation equations for the momentum (1), mass (2), and energy (3), as well as the advection and reaction Equation 4 for the different compositional fields. The energy equation includes the radiogenic heating, shear heating, and adiabatic heating.
Although the model is incompressible, we adopt the equation of state of Murnaghan (Equation 5) (Murnaghan, 1944) to simulate realistic phase transformations that require a temperature and pressure dependent compressible density formulation. Previous studies have shown compressibility to have a small effect on mass conservation for subduction models, suggesting that it can likely be neglected (Fraters, 2014).
where is the final density and refi is the reference density for each composition at surface pressures and a surface temperature of 20°C (Tref), is the thermal expansivity, is the isothermal compressibility, and is the isothermal bulk modulus pressure derivatives.

Rheology
Material in the model has a visco-plastic rheology (Glerum et al., 2018). The viscous regime is handled using a harmonic average of the contribution of dislocation and diffusion creep (Equation 6), whereas the plastic regime uses the Drucker-Prager criterion when the viscous stress exceeds the yield stress (Equation 7). dif f|disl = 0.5 where A is the prefactor rescaled from uniaxial experiment, n is the stress exponent, d and m are the grain size and grain size exponent, e is the square root of deviatoric strain rate, Q is the energy of activation, V is the volume  Oncken et al., 2006), overlain with the extent of the active magmatic arc (red) and the foreland areas with thin-skinned (yellow) and thick-skinned (light-blue) deformation. Blue shaded areas indicate the neighboring flat-slab regions. White arrows show the present day absolute plate velocity . (b) Schematic tectonics of the Altiplano transect at 21°S (dashed rectangle in a), modified from Oncken et al. (2006) and Armijo et al. (2015). The question mark indicates an unclear presence of the lithosphere. (c) Estimated shortening rate evolution (Anderson et al., 2018;Oncken et al., 2006Oncken et al., , 2012, volcanic activity , paleoelevations (Garzione et al., 2017), absolute velocity of the Nazca plate (Sdrolias & Müller, 2006), and perpendicular and parallel convergence velocity (Quiero et al., 2022). Blue squares indicate the margin of error of the estimated retroarc shortening rate (EC + SR, Anderson et al., 2018, modified from Quiero et al., 2022. WC: Western Cordillera; AP: Altiplano Plateau; EC: Eastern Cordillera; SR: Subandean Ranges.
10.1029/2022JB025229 4 of 21 of activation, P the pressure, R the gas constant, and T is the temperature. Here n = 1 for the diffusion creep case and m = 0 for the dislocation creep case. The yield stress is defined by Drucker-Prager law.
where C is the cohesion, P the pressure, and F the internal friction angle in radian. We also included linear plastic strain softening of the friction and cohesion that depends on the strain accumulation over time (Table S1 in Supporting Information S1). For example, in model M1, the internal friction angle of the upper crust decreases linearly from 30° to 6° between an accumulated plastic strain weakening interval of 0-1.5, while the weak foreland sediments decrease from 30° to 3° between 0 and 0.5.
The effective plastic viscosity is calculated by

Numerical Model Set Up
We split the model box into two sub-boxes; a 96 km thick (in depth) box that represents the lithosphere, and an 804 km thick box that represents the sublithosphere mantle. This allows us to set more flexible independent side boundary conditions for each sub-box. For example, during the initialization we prescribe horizontal velocity to the lithosphere, the vertical velocity is stress free, whereas the initial lithostatic pressure is used to simulate an open side boundary in the sublithosphere mantle. The left border is fully open after initialization. The entire box is 2592 × 900 km, which gives an aspect ratio of ∼1:3 while the cells are square (Gerya et al., 2019). The adaptive mesh is refined based on the compositional fields and the magnitude of the strain rate. The asthenosphere mantle and the oceanic lithosphere mantle have a fixed resolution of 32 and 4 km, respectively. Surface topography is calculated using the ASPECT-Fastscape coupling (Bovy, 2021;Braun and Willet, 2013;Neuharth, Brune, Glerum, Heine, & Welford, 2021;Neuharth, Brune, Glerum, Morley, et al., 2021) using a very low diffusion coefficient (∼1e−6 m 2 /year) simulating sluggish surface erosion in the central Andes.
In order to convert the model time (t mod ) to a real time (t real ) we assume that t real (Ma) = 38t mod (My). As such, after the initial flat slab stage (∼6.5 My) the model starting time corresponds to a t real value of 31.5 Ma ( Figure 5). Subduction in the model is initiated by prescribing an oceanic plate velocity of 7 cm/yr in the first 6.5 My, which represents the plate velocity between 35 and 30 Ma (Sdrolias & Müller, 2006). After this no oceanic plate velocity is prescribed, and the oceanic plate freely sinks through the mantle due to slab pull. The trenchward velocity of the continental plate is set to 2 cm/yr, corresponding to the average overriding plate velocity during the last 40 Ma. As gaps in the Andean volcanic activity at ∼30-35 Ma suggest a phase of flat slab subduction (Barazangi & Isacks, 1976;Isacks, 1988;James & Sacks, 1999;Ramos et al., 2002;Ramos & Folguera, 2009), we initialized the model with a flat-subduction stage (Figure 2a, Text S1 in Supporting Information S1). After initialization (Figure 2c), the flat slab segment is ∼250 km long at ∼100 km depth, similar to the current Pampean flat slab (Marot et al., 2014;Rodriguez Piceda et al., 2021). This initialization (see text in Figure S1 in Supporting Information S1) is required in a 2D model to simulate change of buoyancy of the oceanic plate in the Central Andes linked to the southwards migration of the flat slab at ∼25 Ma (Bello-Gonzalez et al., 2018;Yañez et al., 2001). Note that alternative conceptual models also exist. For instance, Kay and Coira (2009)  The geometry of the continental plate is based on structural reconstructions and crustal balance estimations during the Oligocene (Armijo et al., 2015;Hindle et al., 2005;Sobolev et al., 2006). For the shortening analysis, we differentiated two continental domains: the orogenic and the foreland with the thicker lithosphere of the Brazilian Shield margin. We used an oceanic lithospheric thickness consistent with a 40 My old (Maloney et al., 2013) plate near the trench (Turcotte et al., 2002). We assumed a conductive geotherm for the lithosphere and an adiabatic temperature profile for the asthenosphere (Figure 2b) and let the temperature re-equilibrates during the initialization phase.

Generic Models
We have tried to incorporate the most important ingredients found in the literature to simulate the deformation of the Central Andes. In summary, five key ingredients are used to simulate plate interaction in Model 1 (reference model). First, a high-resolution (1 km) visco-plastic subduction interface with a low effective friction coefficient (0.05) enables the brittle-ductile transition to occur at ∼45 km depth. Second, simulate the main phase transitions, Olivine-Wadsleyite-Ringwoodite-Post Spinel transitions for the asthenosphere and lithospheric mantle (Arredondo & Billen, 2016Faccenda & Dal Zilio, 2017) and Gabbro-Eclogite-Stishovite phase transitions for the oceanic crust and continental lower crust to simulate eclogitization (e.g., green color gradient, Figures 4 and 6) and delamination. Third, the rapid weakening of foreland sediments (see Section 2 for details) to allow a transition from a thick-skinned to thin-skinned deformation style, and to initiate underthrusting of the Brazilian cratonic shield Sobolev et al., 2006). Fourth, the prescribed trench-ward motion of the overriding plate velocity, which provides the main driving force for building the Andes (Husson et al., 2012;Martinod et al., 2010;Silver et al., 1998;. Fifth, the flat slab subduction which helps to initiate the thermomechanical weakening of the overriding plate through scrapping of the sublithospheric mantle and its removal, exposing the continental crust to the warmer asthenosphere after steepening of the flat slab (Isacks, 1988;Liu & Currie, 2016). We ran nine alternative simulations to Model M1 (Table S2 in Supporting Information S1): (a) three models with variable interplate friction coefficient ( , Movies S10 and S11). In addition, we provide a model without flat subduction ( Figure S10 in Supporting Information S1, Movie S12, Text S2 in Supporting Information S1) where we illustrate its role; in this context we also discuss the resulting subduction velocity. In all models, we measured the balance between the rate of trench retreat (V tr ), which is positive when the trench migrates westwards, and the overriding plate shortening rate (V sr ), which includes the orogenic shortening rate (V os ) and the rate of underthrusting of the Brazilian cratonic shield (V und ; Table S3 in Supporting Information S1). The orogenic shortening and the underthrusting rate are equivalent to the Interandean and Subandean shortening rate (Oncken et al., 2012). All of these components contribute to accommodating the westward velocity of the overriding plate (V op = V os + V und + V tr ). When the trench retreat rate is less than the overriding plate velocity, the shortening rate increases to maintain the balance. In this context, we refer to the trench as hindered or blocked.

Reference Model (M1)
After the initialization stage, subduction evolved dynamically from 7 to 11 My ( Figure 3a). The slab steepens and accelerates, slowing down the seaward migrating trench as it retreats. Part of the continental mantle starts delaminating when plastic strain localizes over a main thrust fault in the top of the continental crust. During this period, the topographic uplift is limited to the central area or the orogenic domain (equivalent to the Altiplano plateau). At ∼10 My, the block of continental lithosphere consisting of eclogitized lower crust and mantle delaminates and sinks with the slab. At ∼10.5 My, the subduction velocity decreases as trench retreat reinitiates.
From ∼11 to ∼20 My (Figure 3b), relatively fast slab rollback continues as the slab sinks into the transition zone. At ∼18 My, the slab reaches the lower mantle but does not immediately penetrate into it, instead it is deflected  and slowly traverses horizontally along the 660-km phase transition boundary. At ∼20 My, the slab buckles by folding twice to the west and to the east at the transition zone as the trench continues to retreat.
At ∼23.5 My, the slab segment in the upper mantle steepens and halts trench retreat such that the trench no longer moves relative to the mantle. Simultaneously, subduction velocity increases, the previous thrust fault is reactivated, the strain localizes in the eastern orogenic domain, and the lithospheric mantle successively delaminates in the east as the deformation intensifies and migrates east toward the foreland (Figure 3c). Underthrusting of the cratonic shield initiates at ∼26 My during the delamination period. The eastern domain uplifts from ∼20 to 24 My, then slightly subsides at ∼24 Ma.
From ∼25 to ∼31 My the topography significantly uplifts and approximately reaches elevations of the present day (Figures 1c and 6). At ∼29 My, active deformation in the foreland decreases and trench retreat reinitiates as the new slab segment reaches the lower mantle transition trenchward of the older and stalled slab segment. After this time, topography no longer significantly changes ( Figure 5). At ∼30 My, the slab buckles a second time followed by another stage of hindered trench motion at ∼35 My ( Figure 3d) as the slab steepens and accelerates. By ∼33.5 My, the cratonic shield has re-initiates underthrusting beneath the orogenic domain. At ∼37.5 My foreland deformation becomes less efficient and the mantle wedge starts to delaminate as trench retreat reinitiates. Overall, after 38 My the trench retreats ∼330 km, the orogen shortens ∼195 km. Because of underthrusting the foreland shortens by ∼105 km (Table S3 in Supporting Information S1, Figures 9a and 9b, Movie S2).

Models With Variable Interplate Frictions (M2a-c)
We ran three variations of the reference model M1 that has a friction coefficient at the subduction interface of 0.05 (Figures 4b-4d, Movies S3-S5) where we varied the friction: 0.015 (model M2a), 0.035 (model M2b), and 0.06 (model M2c). The friction used in model M2a is thought to be similar to the Southern Andes , although the slab geometry and structure of the upper plate may vary latitudinally. Shortening in model M2a is ∼100 km, with most of the deformation being accommodated within the orogen. No underthrusting occurs in that model, suggesting that the deformation did not reach the foreland. With a friction of 0.035 (model M2b) the orogenic shortening increases (∼170 km), but once the deformation has migrated to the foreland, the shield underthrusts by ∼90 km at ∼29 My. Finally, with a higher friction (0.06) at the interface in model M2c, underthrusting of the shield occurs sooner ∼25 My and reaches ∼105 km. Initially, most of the shortening is  accommodated by the orogen before it quickly migrates to the foreland. The models also suggest that lower friction at the interface results in higher oceanic plate velocities ( Figure S2 in Supporting Information S1).

No Eclogitization Model (M3)
We ran one model (Figure 4e, Figure S7 in Supporting Information S1, Movie S6) without eclogitization of the lower crust to illustrate the importance of this process. As in the reference model, in the model without eclogitization the slab steepens after the flat subduction stage, however, only a small block of the continental mantle is delaminated. Deformation does not efficiently localize as the plastic strain becomes distributed within the orogen. At ∼27.5 My, after the first buckling of the slab, the slab steepens, the orogen shortens, and deformation migrates to the foreland. Shortly after, the shield starts underthrusting as a result of the weakening of the foreland sediments. Underthrusting soon becomes inefficient, but after the second slab buckling (∼34 My), it reinitiates. After ∼38 My, the trench has retreated ∼455 km, the orogenic domain has shortened by only ∼110 km, and the total underthrusting is only ∼65 km.

High Heat Flow Model (M4)
The Central Andes hosts the largest magmatic body in the world (Perkins et al., 2016), and as such the surface heat flux is particularly high (>110 mW/m/K; Hamza et al., 2005;Schilling et al., 2006). This is partly due to partial melting of the crust (up to 20%) as suggested by the high Vp/Vs and high seismic attenuation detected in the area (Haberland et al., 2003;Hamza et al., 2005;Schurr et al., 2003). An increase of the degree of partial melting would result in lower viscosities (Dingwell et al., 1996;McKenzie & Bickle, 1988), and stimulate intra-crustal convection (Arndt et al., 1997;Babeyko et al., 2002). To evaluate the importance of lower viscosities related to partial melt, we ran a model (Figure 4f, Figure S7 in Supporting Information S1, Movie S7) where we increased the thermal conductivity of the upper crust by 1,000× and decreased its minimum viscosity to 2.5e18 Pas (vs. ∼1e22 Pas in the reference model) when the temperature is greater than 1,000 K. In this model the orogenic shortening increases (∼280 vs. ∼195 km in the reference model), with the orogenic domain accommodating the majority of the shortening. Strain in the orogen strongly localizes onto a few faults and, thus, does not migrate to the foreland ( Figure S7 in Supporting Information S1, Movie S7). As a result, the sediments do not accumulate enough plastic strain to weaken and underthrusting does not occur. Shortening rate after the first slab buckling cycle is more efficient than in the reference model and reaches ∼25 mm/yr. Compared to the reference model, the greater orogenic deformation results in a thicker orogen and higher surface heat fluxes (>120 mW/m; Figure S3 in Supporting Information S1).

Foreland Sediments Strength (M5a-b)
Foreland sediments in the reference model (M1) have an internal friction angle of ∼3° and a cohesion of 1 MPa, thus an effective friction coefficient is ∼0.05. We ran two models in which we increased the internal friction angle and cohesion to 10° and 20 MPa (model M5a, Figure 4g, Movie S8) and 30° and 20 MPa (model M5b, Figure 4h, Movie S9). In both models, the model initially evolves like the reference model, but underthrusting is not significant, resulting in a lower accumulated shortening magnitude of ∼240 km ( Figure S8 in Supporting Information S1).

Overriding Plate Velocity (M6a-b)
In the reference model, the overriding plate velocity is 2 cm/yr, which represents an average absolute motion orthogonal to the trench over the last 40 Ma. We ran two alternative models where the overriding plate velocity is 1 cm/yr (model M6a, Figure 4i, Movie S10) and 4 cm/yr (model M6b, Figure 4j, Movie S11). The M6a model results in a total shortening of ∼155 km and a retreat of the trench of ∼160 km ( Figure S9 in Supporting Information S1). The slab piles in the transition zone (Figure 4i); the orogenic domain shortens during the steepening of the flat slab and after its delamination during the last 7 Ma. In model M6b, the slab does not deform and is anchored to the lower mantle, the amplitude of shortening is ∼450 km and the trench retreat is ∼810 km; most of the shortening comes from early foreland underthrusting at ∼23 Ma which stops in the last ∼5 Ma.

Discussion
Our results suggest that the timing of the shortening events is a direct consequence of the interaction between the buckling subducting plate and the weakened overriding plate. In the reference model, we distinguish four notable deformation phases that correspond in amplitude, timing, and space to the shortening rate from the geological compilation (Oncken et al., 2012). Overall, deformation migrates across the orogenic domain to the eastern foreland in four phases, illustrated in Figure 5. The compressive stress generated by the difference of velocity between the trench and the overriding plate is accommodated in one of two ways: (a) orogenic shortening and (b) underthrusting of the foreland. The effectiveness of deformation localization depends on the strength of the overriding plate and the interplate coupling. Here, we discuss the key processes that affect the strength of the overriding plate, the subduction and deformation dynamics of the slab, and, finally, the interaction between the two plates.

Delamination
Extensive lithospheric delamination is known to have taken place under the Altiplano-Puna plateau (Beck & Zandt, 2002;Beck et al., 2015;Kay & Kay, 1993) and contributed to present-day elevations (Garzione et al., 2006(Garzione et al., , 2008(Garzione et al., , 2017Wang et al., 2021). This process is thought to be the result of the eclogitization of the mafic lower crust and lithospheric mantle, which is facilitated by the hydration of the sub-lithosphere from the ∼200 Ma subduction history that accelerates the metamorphic reaction (Babeyko et al., 2002. Additionally,the thick (∼45 km) initial crust at ∼30 Ma results in a high lithostatic pressure in the lowermost crust (Armijo et al., 2015;Hindle et al., 2005;Sobolev et al., 2006). After eclogitization and delamination the crust warms up, which enhances the weakening of the overriding plate and leads to localized deformation and subsequent delamination events. Model M3 demonstrates that delamination and shortening are inhibited without eclogitization (Figure 4). Whereas model M4 shows that, due to the faster crustal thickening accumulated by the weak crust, eclogitization and delamination are very effective when the orogenic domain is thermally weakened. In that latter case, the migration of the deformation to the foreland is not guarantee ( Figure S7 in Supporting Information S1).
We observe two delamination stages after the first event caused by steepening of the flat slab in Phase I. First, the initial removal exposes the crust at the western edge that is directly in contact with the asthenosphere, thereby increasing its temperature and decreasing the viscosity at its base. As a result, the lower crust delaminates faster in the west, causing it to asymmetrically drip to the east (i.e., Stage 1 in Figure 6a). The pure shear deformation localizes in the orogenic domain until delamination is complete. Second, when the viscous deformation of the orogen connects with the plastic deformation of its foreland at 26 My, the foreland start underthrusting beneath the orogen due to the low effective friction of sediments. This results in orogenic thickening and a switch from pure shear to simple shear shortening. Consequently, deformation migrates to the east causing delamination to accelerate (Stage 2 in Figure 6b). We note that our model contrasts with the previous model of Sobolev et al. (2006) where the delaminating lithospheric mantle flowed toward the subduction wedge, and coupled with the sinking plate, increasing the shortening rate.

Mechanical Weakening of the Foreland Sediments
The presence of weak foreland sediments in our model is the key factor in simulating the transition from pure shear deformation to simple shear deformation at ∼10 Ma in the Altiplano. Simple-shear shortening is associated with higher strain localization over fewer faults and the formation of deep low-angle detachments. In the foreland, these faults are situated near the base of the sedimentary cover and are characteristic of the thin-skinned tectonic style. Increased fluid pressure in the Paleozoic shale layers (Allmendinger & Gubbels, 1996), likely due to rapid deposition of foreland basin strata (Uba et al., 2009), at the front of the orogen, may have resulted in transient weakening and reduction of the effective coefficient of friction to ∼0.05 or less, initiating the underthrusting of the Brazilian cratonic shield Babeyko et al., 2006). In the reference model, underthrusting takes place in two stages. The first stage happens during hindered trench motion at ∼20.5 My, causing the deformation to migrate to the foreland. When the active brittle shear zone, from the failure of the foreland sediments, connects to the ductile shear zone accommodating the on-going delamination, underthrusting becomes more efficient. The delamination also facilitates the underthrusting of the Brazilian cratonic shield that meets less resistive forces. Underthrusting of the shield forces the middle and lower crust to flow and thicken forcing the topography to uplift, reaching present-day elevations of ∼4 km at ∼31 My (∼7 Ma ago). A second stage of underthrusting occurs in the last ∼4 My when the trench is again blocked, but the topography does not change significantly ( Figure 5).

Subducting Plate
While the westward motion of the South American plate provides the main force (Husson et al., 2012;Martinod et al., 2010) for the tectonic shortening, the magnitude of the compressive stress in the South American plate margin is determined by the resistance of the Nazca plate (i.e., by the ability of the trench to retreat; Funiciello et al., 2008;Holt et al., 2015;Lallemand et al., 2005Lallemand et al., , 2008. In the Central Andes, the trench has migrated west over the last ∼40 Ma as a result of the rollback and subsequent sinking of the bending slab in the asthenosphere, as well as the forced trench retreat from the excess velocity of the overriding plate (Schepers et al., 2017). Recent studies have proposed that the trench velocity can also be affected by deep subduction dynamics (Boutoux et al., 2021;Briaud et al., 2020;Faccenna et al., 2017). In this section, we discuss the implications of these subduction dynamics.

Flat Slab Steepening
The cause of flat subduction is still debated. It likely results from the shallowing of the slab from long lasting subduction, as well as larger buoyancy related to the Juan Fernandez ridge (Schellart, 2020;Schellart & Strak, 2021) that has migrated to the south during the last ∼35 My (Figure 1; Bello-González et al., 2018;Yáñez et al., 2001). Others authors (Liu & Currie, 2016;Quinteros & Sobolev, 2013) suggested that slab break-off in the upper mantle could have further contributed to the slab flattening by decreasing the slab pull force. The flattening of the slab could also be caused by an increase in asthenospheric pressure due to the proximity of a thick cratonic lithosphere (Manea et al., 2012(Manea et al., , 2017Pérez-Gussinyé et al., 2008).Most shortening in the Central Andes occurs after the passage of the ridge (Oncken et al., , 2012, so in this study we focus on processes after the flattening event. Our models suggest that a flat slab at ∼100 km depth, analogous to the Pampean flat slab, could scrape the base of the lithosphere. Eventually at ∼7 My model time, the slab steepens and accelerates as the trench becomes blocked (Figure 7a). The continental mantle coupled to the flat slab segment blocks the corner and is pulled down and thus accommodates the deformation. When the lower continental crust eclogitizes, plastic strain localizes in the top portion of the crust. Slab steepening then accelerates due to the eclogitization of the oceanic crust and parts of the lithosphere are removed. This process of delamination is similar to the mechanism of "blocking of the subduction corner" of Sobolev et al. (2006) for which the increase in shortening rate results from the coupling between the delaminated lithosphere driven downward by the slab. According to our alternative model M7 ( Figure S10 in Supporting Information S1) that has no flat slab, flat slab steepening plays a key role in triggering the initial weakening of the overriding plate, and is facilitated by lower-crustal eclogitization.

Buckling Instability Cycles
Slab buckling occurs when the oceanic plate subducts into the more viscous mantle transition zone or the lower mantle. The difference in velocity between the deeper slab segment relative to the new subducting segment, is accommodated by slab deformation (Ribe et al., 2007). Previous studies have suggested that the lower mantle viscosity and the dip, age, thickness, and strength of the oceanic plate may affect the buckling periodicity and timing of slab stagnation in the transition zone, and additionally could be linked to periodic crustal deformation (Boutoux et al., 2021;Briaud et al., 2020;Capitanio et al., 2011Capitanio et al., , 2010Cerpa et al., 2014;Čížková & Bina, 2013;Garel et al., 2014;Lee & King, 2011;Lyu et al., 2019;Marquardt & Miyagi, 2015;Quinteros & Sobolev, 2013;Quinteros et al., 2010;Ribe et al., 2007). Analyzing the variety of interchangeable parameters that affect the buckling process exceeds the scope of this study. Here, we first interpret the different stages of the buckling cycles and then propose that the westward velocity of the upper plate is a primary factor in controlling the subduction dynamics regime.
We identified two buckling cycles, at ∼20 My and at ∼30 My. Within each cycle, three main events are distinguished that may affect the trench migration rate: 1. Slab impediment (Figure 7b) takes place when the slab meets viscous resistance. This is the case when the slab is impeded by the viscous lower mantle at the beginning of a buckling cycle (∼17 and ∼29.5 My), or before steepening. For instance, when the slab reaches the viscous lower mantle it does not immediately penetrate it. The first slab segment in contact with the lower mantle slows down and viscously resists the new, still sinking, segment. This difference of velocity between the two segments is accommodated through bending in the slab. During these slab impediment events the dip of the slab becomes shallower and the trench continues retreating. This mechanism differs from slab anchoring (Faccenna et al., 2017), in which the difference of velocity between the two segments is too small to cause the folding of the slab. 2. Slab folding (Figure 7c) events occur when, after slab impediment, the slab dip flips in the transition zone.
The now shallower slab dip enables the trench retreat, though no significant deformation is observed. Each buckling cycle consists of two folding events, each consisting of a syncline and an anticline at ∼20, 21 My and ∼30, 33 My, respectively. 3. Slab steepening (Figure 7d) is a drastic event that occurs at the end of a buckling cycle after the second folding event (∼23.5 and ∼33.5 My). Chronologically, the sinking slab meets resistance from the last fold to the east (i.e., Impediment) and bends to the west as for the first folding event. However, the overriding plate has forced the trench to retreat during the previous events, which, prevents the slab from piling up. The slab continues to sink in the transition zone, steepens, and accelerates. The trench slows down and blocks the overriding plate that shortens to accommodate the ongoing deformation. When the trench is blocked, the horizontal stress in the overriding plate can reach values of ∼350 MPa ( Figure S1 in Supporting Information S1, movie S1), which exceeds the maximum strength of the crust (∼250 MPa, Figure 2) and causes it to shorten. Overall, slab shallowing is associated to periods of trench retreat related to the folding events, whereas slab steepening is associated to periods of hindered trench motion following folding events.
This chain of events occurs in a single subduction dynamics regime, primarily determined by the absolute motion of the overriding plate orthogonal to the trench. By comparing the subduction velocity to the trench velocity, we can identify three regimes of subduction dynamics. In the first regime, the lower velocity of the overriding plate leads to the "piling" of the slab in the transition zone (Regime 1; Model M6a, Figures 4i and 8b and Figure S9 in Supporting Information S1). In the second regime, the trench episodically blocks due to "buckling and steepening" of the slab (Regime 2; Model M1, Figures 4a and 8a and Figure S10 in Supporting Information S1). In the third regime, the "anchoring" of the slab in the lower mantle and the high forced trench retreat prevent it from buckling (Regime 3; Model M6b, Figures 4j and 8c and Figure S9 in Supporting Information S1).
In all models, flat slab steepening occurs in a similar manner (i.e., Phase 1). However, westward migration of the trench (i.e., Phase 2) is associated with steady trench retreat for M1 and M6a, but not for M6b for which delamination began earlier. In the latter case, the velocity of the overriding plate is absorbed by the deformation of the orogenic domain, its early delamination, and the underthrusting of the foreland, as the craton no longer encounters resistance from the lithospheric mantle. The rates of shortening and of trench retreat further intensify when the slab becomes anchored in the mantle and the craton is blocked by the slab (Figure 8c and Figure S9 in Supporting Information S1). Both models M1 and M6a are characterized by slab folding, but the higher trench  Lallemand et al. (2008) for North Chile (NC) and Peru (P) using different reference frame NNR-NUVEL1A and HS3-NUVEL1A (NNR and HS3, Gripp & Gordon, 2002) and SB04 (Steinberger et al., 2004). The gray line is the regression  for which subduction zones at present day trend to align with using HS3. retreat of M1 allows episodic slab steepening, resulting in trench blocking and a slight acceleration of the subduction rate. In contrast, in model M7 that we ran without a flat slab ( Figure S10 in Supporting Information S1), no deformation happens implying that buckling and steepening occur independently. Nevertheless, the initial weakening of the overriding plate caused by the passage of the flat slab is necessary to trigger the shortening. Thus, the strength of the orogenic domain is controlled by the timing and intensity of delamination, which plays a major role in transmitting the velocity of the overriding plate to the trench, and controls the timing of underthrusting. The trench velocity ultimately determines the regime of subduction dynamics. Given that westward movement of the overriding plate has decreased from 45 Ma (∼3 cm/yr) to the present day (∼1 cm/yr), we suggest that the Andean subduction regime may have changed from anchoring to buckling and steepening over the last ∼20 Ma.

Interplate Coupling
Our models predict that an effective friction of 0.035-0.05 is required in the Central Andes to obtain significant deformation that is consistent with previous estimates Sobolev et al., 2006). Higher friction values result in slightly slower subduction velocity ( Figure S2 in Supporting Information S1) but more intense pulsatile shortening phases during slab steepening ( Figure S6 in Supporting Information S1). The effective friction is dependent on the sediment thickness at the trench, which at present day may vary from ∼0.5 to ∼2 km in the Central and Southern Andes, respectively (Lamb & Davis, 2003). This latitudinal variation results from the efficiency at which the surface processes supply sediments to the trench. In the last ∼6 Ma, glacial erosion supplied a large amount of sediments to Southern Andes trench. Whereas in the Central Andes, the internal drainage of the Altiplano-Puna plateau is related to low erosional rates that have contributed to sediment starvation at the trench (Lamb & Davis, 2003;Melnick & Echtler, 2006;Oncken et al., 2006). Hu et al. (2021) also show that the temporal variation of sediments included in subduction can lead to a temporal variation in the interplate coupling, suggesting that sediments can be a major contribution of the along-strike shortening magnitude and, thus, the onset of deformation.

Slab Buckling and Overriding Plate Interaction
The rapid growth of the Andes in the last ∼10 Ma (Figures 1c, Garzione et al., 2017) results from a sequence of events generated by plate interactions. While subduction dynamics exerts a major control on plunging plate deformation by blocking trench migration, the overriding plate strength ultimately controls where deformation localizes and forces the trench to retreat. In model M1, when the slab does not steepen, the trench is forced to retreat at the prescribed westwards velocity of the overriding plate, ∼2 cm/yr (M1, Figure 9c). Alternatively, in model M3 without lower crustal eclogitization, the slab steepens independently of the shortening of the stronger overriding plate and the velocity of the overriding plate is accommodated by the forced trench retreat ( Figure S7 in Supporting Information S1). This latter case indicates that the steepening is mostly controlled by slab strength and the slab buoyancy rather than the shortening of the overriding plate. The upper plate strength is evolving, first, with the passage of the flat slab that may have initially weakened the lithosphere through partial removal of the lithospheric mantle, and through thermal weakening related to crustal exposure near the hotter asthenosphere (Isacks, 1988), and second, by triggering the subsequent delamination (see previous section).
Pulsatile behavior in the deformation of the South American plate has been inferred from paleoelevation reconstructions using stable isotope (Boschman, 2021;Garzione et al., 2008;Leier et al., 2013) and the magmatic activity (Decelles et al., 2009). We suggest that buckling instabilities in a subducting plate offer a plausible explanation in the variability and timing of the Nazca plate deformation during the last ∼20 Ma as well as the present-day deep seismicity distribution ( Figure S4 in Supporting Information S1). We find that shortening rate pulses occur at the end of each buckling cycle when slab steepening inhibits trench retreat (Figures 9c and 9d), and that these pulses reproduce similar signals to what is seen in the geological data. When the slab steepens the forced trench retreat from the overriding plate is hindered and the horizontal stresses increase ( Figure S1 in Supporting Information S1; Movie S1), resulting in a shortening of the upper plate. Additionally, in the last ∼2 Ma the geological data shows a decrease in the shortening rate, which is also predicted by our model through underthrusting. At later stages, the trench retreat resumes and underthrusting loses its efficiency, which could indicate the beginning of a new buckling cycle. More recently, based on updated high-resolution convergence rate data orthogonal to the trench (Figures 1c), Quiero et al. (2002) show that there are some short-term variations over the last ∼30 Ma.
They attribute these variations to the delamination of the Central Andes (Quiero et al., 2022). The model M7 with no flat slab suggests that these variations could represent pulses associated with deep subduction dynamics, causing a periodic increase in subduction velocity ( Figure S10 and Text S2 in Supporting Information S1). Nonetheless, in this model, no deformation occurs due to the lack of weakening of the overriding plate from the absence of the flat slab and the overriding plate remains strong, so trench retreat is therefore more effective and may delay the steepening of the slab, causing the last pulse to occur ∼3 Ma later ( Figure S10d in Supporting Information S1). We also notice that the periods of steepening also correlate with flare-up of volcanic activity and greater volume of ignimbrites ( Figure S1c in Supporting Information S1, Trumbull et al., 2006). Slab steepening and lithospheric mantle delamination becomes more active when the trench is hindered, which can lead to ( Figure 9a) an increase in the basal heat flow of the lithosphere and more intense volcanic activity (Section 4.1.1, Isacks, 1988;Kay & Coira, 2009).
Previous seismic tomography studies indicate two large negative seismic anomalies near the transition zone (at depths of 600 and 900 km) that are attributed to slab accumulations (Chen et al., 2019;Liu et al., 2003;Widiyantoro, 1997). The deeper accumulation may relate to a slab anchoring before (Faccenna et al., 2017, Figure  S4 in Supporting Information S1), suggesting that previous accumulation cycles could have occurred before and have "avalanched" in the lower mantle (Briaud et al., 2020;Hu & Gurnis, 2020), wherein they may have become detached from the shallower slab. Indeed, the Peruvian phase from ∼110 to 70 Ma consisting of a rapid series of alternating compressive and extensive regimes of ∼10 Ma (Faccenna et al., 2017;Mora et al., 2009) may indicate that slab buckling events have happened earlier in the northern subduction history. However, because of the absence of an efficient weakening mechanism to trigger delamination and too thin crust to facilitate eclogitization, the orogen experienced no significant deformation. Potentially, we suggest that these avalanche events may have repeated at least three times over the last ∼90 Ma, as suggested by the three cycles of convergence rate recognized in Martinod et al. (2010). These events could also explain the cyclicity of the orogenic processes (Decelles et al., 2009(Decelles et al., , 2014Haschke et al., 2002).

Model Limitations
The main limitation of our model is its dimensionality. The use of 2D modeling is appropriate for the Central Andes, where toroidal flow affecting the edges of the Nazca plate can be neglected. However, latitudinal crustal flow is estimated to have contributed between ∼10% and 30% of the present day crustal thickness of the Central Andes (Hindle et al., 2005;Kley & Monaldi, 1998). In our models the crustal thickening is mainly caused by intraplate shortening. As a result, the crustal thickness of the orogen in our models is lower than the actual thickness of the Central Andes. For example, in model M1 the final orogenic crustal thickness is ∼57 km, whereas it should increase to ∼62-74 km (+10-30%) taking into account the latitudinal crustal flow. In addition, shortening may also be underestimated, which could partially explain the crustal thickness deficit in the models.
In the model M2a ( Figure S6 in Supporting Information S1), the final dip of the slab agrees with the seismic tomography ( Figure S4 in Supporting Information S1). In this model the interplate friction is 0.015, similar to the expected friction of the southern Andes . This suggests that the current dip of the subducting plate in the central Andes is partially caused by lateral support of the shallower oceanic plate to the south. In addition, the appearance of a deeper slab pile 900 km to the east could also indicate deep mantle flow that is not accounted for in our model. Slab buckling could provide a plausible explanation for the distribution of deep seismicity ( Figure S4 in Supporting Information S1).
We find that with an interplate effective friction of 0.05 ( Figure S2 in Supporting Information S1), the maximum amplitude of the modeled subduction velocity is lower than the absolute normal velocity of the Nazca plate (∼12.5 cm/yr at ∼20 Ma, Sdrolias & Müller, 2006). This suggests that slab velocity is largely controlled by the slab interface friction southward and northward of central Andes, where effective friction is lower due to the presence of sediments. Interestingly, despite the lower subduction velocity, our model predicts a decrease in the pulsatile intensity during the last ∼25 Ma to 20 Ma in accord with the observations.

Conclusion
In this study, we demonstrate that cycles of slab buckling due to dynamic slab behavior can explain the timing and amplitude of tectonic shortening pulses observed in the Central Andes since the Late Eocene. The findings of our subduction-related Andean models, suggest that the primary cause of these pulses that contributed to the growth of the Central Andes is the evolving geometry of the subducting Nazca plate. In particular, the steepening of the slab in the upper mantle slows down the trench retreat and subsequent shortening of the advancing South American plate. This steepening first occurs after the end of the flat slab episode at ∼25 Ma. By eroding the lower part of the mantle lithosphere, this episode predisposes the margin for the next deformation phases by decreasing its strength. Later, slab steepening occurs following the buckling of the slab in the mantle transition zone. This new buckling-steepening mechanism sheds light on the causes of the rapid pulsatile growth of the Central Andes during the last ∼20 Ma, and the model evolution is consistent with geological data (Oncken et al., 2012) and with the timing of uplift of the Altiplano plateau (Garzione et al., 2017). Our study also confirms the previous modeling results Liu et al., 2022;Sobolev et al., 2006) regarding the important roles of long-term overriding by South America plate, high intraplate friction due to the lack of sediments in subduction channel, lithosphere delamination of the lithospheric mantle, and weakening of the foreland sediments in the shortening evolution of Andes.