Plateau Formation Controlled by Lithospheric Foundering Under a Weak Crust

Lithospheric foundering is hypothesized to contribute to the formation of Earth's largest continental plateaus, but the predicted effects remain poorly constrained, especially considering variations in crustal strength. Here we propose that lithospheric foundering involving a hot, weak crust can explain aspects of topography and crustal deformation in mountain belts. We explore this hypothesis using numerical models of orogenesis and foundering in the Southern Puna Plateau in the Central Andes. Contrary to previous models of foundering involving a strong crust, which are characterized by subsidence and shortening, a weak crust results in surface uplift and upper‐crustal horizontal extension, which is accommodated by horizontal shortening in peripheral regions. Our model explains data such as the timing and location of exhumation and basin sedimentation as a response to foundering. The model also shows that foundering is capable of converting a high‐relief, broken‐foreland region into a high‐elevation, low‐relief plateau in a convergent tectonic setting.

an orogen generated by crustal shortening, thickening, and magmatism due to oceanic plate subduction (e.g., Allmendinger et al., 1997). In the southern part of the Central Andes, the Puna Plateau ( Figure 1) hosts a thick crust of 55-65 km, but seismic and gravimetric evidence suggests that the mantle lithosphere is only 5-10 km thick at present, suggesting it has been almost completely removed (Heit et al., 2014;Tassara et al., 2006). The Central Andean crust has been deformed multiple times since the Paleozoic, and was greater than 45 km thick at the start of the Cenozoic orogenesis (Profeta et al., 2016). The Southern Puna Plateau (SPP) also hosts a Late Miocene-Recent volcanic province with monogenetic mafic lavas potentially derived from partial melting of foundering lithosphere (R. W. Kay & Kay, 1993;Murray et al., 2015). Elevated crustal temperatures and crustal melting are indicated by silicic volcanism and low shear-wave speeds (Ward et al., 2013). A felsic composition, elevated temperatures, and inherited structures have all contributed to a rheologically weak crust in the Puna (Babeyko et al., 2002).

Results
We present a thermal-mechanical model of lithospheric and asthenospheric deformation ( Figure 2). The model includes kinematic velocity boundary conditions, topographic deformation, and undergoes nonlinear viscoplastic deformation in both the lithosphere and the asthenosphere. We use the open-source geodynamic modeling software ASPECT (Advanced Solver for Problems in Earth's ConvecTion; Gassmöller et al., 2018;Gouiza & Naliboff, 2021;Heister et al., 2017;Naliboff et al., 2020).
The model includes a granite and felsic granulite crust that undergoes tectonic shortening, crustal thickening, eclogitization, and lithospheric foundering. A viscoplastic rheology with strain softening (Table S1 in Supporting Information S1) acts to weaken the crust during deformation. The central portion of the model, which we refer to as the orogen, includes a broad region of damaged upper crust and pre-imposed crustal shear zones ( Figure 2). The model assumes that eclogitization is limited to the central portion of the orogen (Figure 2) such that the width of the eclogitic root matches the ∼65 km width of a positive tomographic anomaly previously identified  (Schnurr et al., 2006) and extensional faults, after Schoenbohm and Carrapa (2015). Moho contours given in km (Tassara et al., 2006). as a foundered lithospheric block (Bianchi et al., 2013;Calixto et al., 2013). We focus our discussion on a base model (v0), but based on our analysis of a suite of models varying key parameters (v1.1-v10.3), all models result in a broadly similar sequence of events ( Figure S3 in Supporting Information S1), which we classify into five distinct temporal phases ( Figure S4 in Supporting Information S1). See Text S1 in Supporting Information S1 for a full model description.

Model Evolution
Phase 1. The first phase (0-22 Myr in model time) is characterized by horizontal shortening of the orogen. Shortening is accommodated by deformation on distinct structures in the upper crust and by viscous thickening of the middle to lower crust and mantle lithosphere. By the end of phase 1, the orogen is characterized by a rugged contractional basin-and-range topography with a mean elevation of ∼2 km and topographic relief of ∼2 km ( Figure 3, 15 Myr; Figure S5 in Supporting Information S1).
Phase 2. The second phase (22)(23)(24)(25)(26)(27)(28)(29)(30) is characterized by the formation of a negatively buoyant root. Isostatic subsidence of the root results in buckling of the lower crust and mantle lithosphere, producing convergent Poiseuille flow and shortening of the middle-upper crust above the root (Figure 3, 30 Myr). This accelerated shortening in the central sector of the orogen outpaces the imposed tectonic shortening rate, so a drastic reduction in Figure 2. Numerical approach. (a) Strong lithospheric blocks (shaded gray) converge at a rate v a toward the central "orogen" with initial plastic damage in upper crust ("initial damage," yellow to brown color scale). Small outflux velocities are prescribed to conserve mass (v b and v c ). W ec denotes the width of eclogite-prone lower crust. (b) Eclogitization as a function of temperature and pressure. (c) Temperature (red) and strength profiles at the end of Phase 2, assuming a strain rate of 10 −15 s −1 . Dashed lines are plastic yield strengths for 2° and 15° effective friction angles (DG = dry granite, FG = felsic granulite, DO = dry olivine, Ec = eclogite). Material properties are given in Supporting Information S1 (Table S1).

of 12
shortening rate occurs in the peripheral regions ( Figure 3, 30 Myr); in models with weaker crusts (v3.1) horizontal extension occurs in these regions. The formation of the eclogitic root therefore reconfigures the distribution of deformation within the orogen for approximately 8 Myr, after which the root begins to founder. Phase 3. The third phase (30-37.6 Myr) is characterized by foundering of the eclogitic root as a viscous drip. The transition from isostatic sinking to dynamic growth of the drip is marked by an inflection point in the maximum depth of eclogite ( Figure S4 in Supporting Information S1). The drip induces flow within the crust, which drives both crustal thickening and upper-crustal extension (Figure 3, 37 Myr). Extension is confined to a 50-75 km-wide region centered on the drip and reactivates inherited contractional structures. This horizontal extension is accommodated by accelerated shortening on the margins of the orogen, where shortening had previously been reduced during Phase 2 ( Figure 3, 37 Myr). The mean elevation and topographic relief within the orogen increase during the growth of the drip ( Figure S5 in Supporting Information S1). Surface uplift is greatest in the center of the orogen, which resembles a 75 km wide proto-plateau by the end of Phase 3 ( Figure 3, 37 Myr).
Phase 4. The fourth phase (37.6-37.8 Myr) is characterized by the detachment of the drip, a rapid process that lasts 200 Kyr ( Figure S4 in Supporting Information S1). The crustal welt that was supported by the drip collapses gravitationally, yielding a strong, transient pulse of horizontal extension that initially encompasses most of the orogen (Figure 3, 37.7 Myr). Extension and crustal heating are greatest directly above the detaching lithosphere ( Figure 3, 37.7 Myr). As a result, mountain range elevations and topographic relief rapidly decrease in the center of the orogen ( Figure S5 in Supporting Information S1, and compare Figure 3 at 37.7 and 45 Myr). Despite this extensional collapse, foundering ultimately yields ∼500 m of isostatic uplift of the entire orogen ( Figure S5 in Supporting Information S1). By the end of Phase 4, the central 250 km of the orogen resembles a high-elevation low-relief plateau (Figure 3, 37.7 Myr).
Phase 5. The fifth phase (37.8-50 Myr) is characterized by the gradual return of tectonic shortening to the orogen, referred to here as "recovery." Because of the model setup, the descending drip migrates west relative to the overriding lithosphere. The westward migration of the drip drives upwelling flow in the asthenosphere that preferentially thins and heats the lithosphere east of the center of the drip ( Figure S6 in Supporting Information S1). A horizontal offset of 20-100 km (50 km in the base model) develops between the drip and the maximum in crustal temperature ( Figure S6 in Supporting Information S1). Previous studies show that mantle corner flow is capable of pulling a descending drip trenchward (i.e., westward) Wang et al., 2021), and the 50 km offset in our base model is similar to such results (Wang et al., 2021).
The asymmetry in drip detachment produces heterogeneities in crustal thickness and temperature, contrasts in gravitational potential energy, and differing rates of recovery for the western and eastern half of the orogen. In the western half, crustal shortening resumes ∼4-5 Myr after drip detachment (Figure 3, 45 Myr). In the eastern half, thermally buoyant crust resists shortening and is instead placed into extensional to neutral conditions for ∼7-12 Myr (Figure 3, 45 Myr). This extensional-neutral region reduces in width at a rate proportional to tectonic convergence rate until the entire orogen experiences horizontal shortening by ∼50 Myr (Figure 3, 50 Myr). This kinematic reorganization results in a nearly uniform crustal thickness of 60 km and a mean elevation of 4.3 km ( Figure 3, 50 Myr; Figure S5 in Supporting Information S1).

Discussion
Although our models show promising results that can be applied to natural orogens, care is required to interpret them in the context of the Central Andes. Lack of deformation in the 2-D models could be consistent with out-of-plane deformation (strike-slip) in a 3-D context. Our models parameterize the effects of mantle corner flow on a lithospheric drip, but they do not include an oceanic slab, which would strongly affect the flow and viscosity of the asthenosphere. The 100-150 km deep Nazca slab beneath the Southern Puna could have interacted with foundering lithosphere, as suggested by some tomography studies (Calixto et al., 2013). These and other important limitations are discussed in Supporting Information S1 (Text S1).

Comparison to the Southern Puna
In the SPP, lithospheric foundering was originally proposed based on the co-occurrence of horizontal extension and back-arc mafic volcanism (R. W. Kay & Kay, 1993;S. M. Kay et al., 1994). Extension and volcanism both initiated during the Middle to Late Miocene in the Antofalla region  (the central sector of the SPP) and migrated outward toward the edges of the plateau during the Pliocene (Schoenbohm & Carrapa, 2015;Schoenbohm & Strecker, 2009). The volcanics carry an isotopic signature of garnet-pyroxenite melting, suggesting the presence of an eclogitic crustal root prior to ∼7 Ma (Ducea et al., 2013;Murray et al., 2015).
Sedimentological and structural studies also place the detachment of the SPP lithosphere ∼7 Ma (Kraemer et al., 1999;Schoenbohm & Carrapa, 2015). The following discussion therefore assumes that drip detachment at ∼38 Myr in model time equates to 7 Ma on the SPP timescale, allowing us to directly compare the model evolution with geologic data (Figure 4).
Existing thermochronometry data from across the plateau suggest a cessation of exhumation in the western (Quebrada Honda Range, McMillan et al., 2022) and eastern margins (Laguna Blanca and Chango Real Ranges) between ∼20 and 12 Ma (22-30 Myr in the model), but these data indicate that shortening continued in the central sector (Calalaste Range and Salina del Fraile) during this time (Carrapa et al., 2005;Kraemer et al., 1999;McMillan et al., 2022;Zhou et al., 2017). Such an evolution of deformation is consistent with the narrowing of the actively deforming area in the model due to eclogitization (Figure 4a).
Horizontal extension initiated between 16 and 11 Ma in the Antofalla region, west of the Calalaste Range, by reactivating a major reverse fault (the Antofalla Fault) as a normal or oblique normal fault Tye et al., 2022). In contrast, horizontal shortening is documented in the surrounding regions of the SPP during this time period (Figure 4c), and outward propagation of shortening both westward and eastward from the Antofalla region is documented. Westward propagation is suggested by thrusting ∼18 Ma on the margin of the Salar de Antofalla (Kraemer et al., 1999) and ∼16 Ma in the Salina del Fraile . Eastward propagation of shortening during this time is well known (Zhou et al., 2017), and is also suggested by the ∼15-10 Ma exhumation of the Laguna Blanca Range  and ∼11 Ma syn-contractional strata in the Pasto Ventura region (Schoenbohm & Carrapa, 2015). The localized initiation of extension and a coeval broadening of the area undergoing shortening shown by the model during Phase 3 therefore seems to be reflected in existing data for the SPP (Figure 4c).
Horizontal shortening resumed in the Antofalla region approximately 3.6 and 2 Ma, as bracketed by normal faulting of 3.6 Ma volcanic flows and overthrusting of 2-3 Ma Quaternary deposits (Adelmann, 2001). Ongoing horizontal shortening beneath the Calalaste Range indicated by upper-crustal earthquake focal mechanisms, which show thrust-sense displacement (Mulcahy et al., 2014). In the eastern sector of the SPP, however, Late Miocene shortening was followed by Pliocene to Recent horizontal extension (e.g., Allmendinger et al., 1989;Marrett et al., 1994;Schoenbohm & Strecker, 2009;. The SPP is therefore characterized by horizontal shortening of the central sector and extensional and strike-slip deformation of the eastern sector ( Figure 4d). As shown by the model, crustal buoyancy contrasts following asymmetric detachment of a lithospheric drip offer a viable explanation for the complex Pliocene to Recent deformation of the plateau. Although our model does not include melting, its crustal heating pattern ( Figure S6 in Supporting Information S1) is consistent with the observed patterns of melting and silicic volcanism on the SPP, which are also reproduced in models of foundering with mantle corner flow (Wang et al., 2021). The model predicts that crustal heating and extensional deformation are co-located along the southeastern margin of the plateau because they are both consequences of asymmetric drip detachment.

Implications
Our results confirm previous models of upper crust extension during drip growth (Conrad & Molnar, 1997;Molnar, 2015) and reveal additional behaviors not previously investigated. Low to moderate extension and surface uplift rates during Phase 3 are consistent with an overall weak crust with a temperature-dependent viscosity that decays with depth (Molnar, 2015). Additional behaviors include extensional collapse of a crustal welt following drip detachment (even within a convergent tectonic setting) (Figure 3 and Figure S5 in Supporting Information S1, 37.7 Myr), differential crustal heating and buoyancy following drip detachment ( Figure S6 in Supporting Information S1), and the development of broad low-relief topography due to extensional reactivation of range-bounding structures (Figure 3 and Figure S5 in Supporting Information S1, 45 Myr).  . Major thrust faults (Schnurr et al., 2006) and normal/strike-slip faults (Schoenbohm & Carrapa, 2015;Tye et al., 2022) are only drawn according to inferred shortening and extending regions, respectively. Apatite fission track (AFT) data (Carrapa et al., 2005;McMillan et al., 2022;Stalder et al., 2020) and syntectonic sediments (Adelmann, 2001;DeCelles, Carrapa, et al., 2015;McMillan et al., 2022;Schoenbohm & Carrapa, 2015; are mapped by age. Panel d shows earthquake focal mechanisms (Mulcahy et al., 2014) and a P-wave anomaly from <100 to ∼200 km depth (Beck et al., 2015). Red line denotes modeled cross-section and major structures labeled in (a).
Our model sheds light on the geodynamic evolution of orogenic plateaus, especially those for which lithospheric foundering has been hypothesized. Specifically, we suggest that foundering-induced crustal flow can drive lower crustal thickening and heterogeneous upper-crustal extension and shortening during orogenesis, a robust outcome predicated largely on the presence of a relatively weak, felsic lower crust (Table S2 in Supporting Information S1). Although the lower continental crust is sometimes assumed mafic, recent analyses strongly question this assumption, especially for major orogens (Hacker et al., 2015;Sammon et al., 2022). We therefore suggest that lithospheric foundering in regions such as the Altiplano (Central Andes) (e.g., Sundell et al., 2019), Tibetan Plateau (e.g., Kapp & DeCelles, 2019), and Western U.S. (e.g., Smith et al., 2017) may have involved significant crustal flow and major episodes of upper crustal extension.
The orogenic processes operating in the Central Andes remain an ongoing topic of research and debate (DeCelles, Zandt, et al., 2015;Gianni et al., 2020;Zapata et al., 2020). It is well known that upper-crustal shortening is not sufficient to account for the 60 km crustal thickness in the SPP (Kley & Monaldi, 1998). Our model suggests that middle and lower crustal flow induced by foundering played a significant role. Although other mechanisms such as lithospheric delamination (Beck et al., 2015) and flat-slab subduction (Ramos & Folguera, 2009) may offer alternative explanations for out-of-sequence deformation and volcanism, it is unclear how they would produce the early onset of extension in the Antofalla region compared to the eastern sector of the SPP . Pliocene-Quaternary extensional deformation of the SPP is often ascribed to orogen-parallel stretching (Allmendinger et al., 1989;Marrett et al., 1994;Riller & Oncken, 2003) or gravitational spreading (Daxberger & Riller, 2015;Schoenbohm & Strecker, 2009), but these processes also do not explain the complex deformation history of the interior of the SPP  and are not mutually exclusive with the foundering hypothesis.
The emplacement of a dense lithospheric root reduces the mean elevation of the orogen and concentrates shortening within the orogenic hinterland. The root's foundering rapidly increases mean elevation and reduces topographic relief within the orogen ( Figure S5 in Supporting Information S1), and facilitates propagation of shortening toward outward toward the foreland (Figure 4). In NW Argentina, both the SPP and the broken-foreland region are deformed by mountain ranges that share a similar short-wavelength topography, the plateau is more elevated, more highly deformed, and has a lower topographic relief . We hypothesize that lithospheric dripping may play a role in setting these topographic differences between the plateau and foreland regions (Figure 1). Foundering may also explain why upper crustal shortening propagated eastward from the central portion of the plateau into the foreland region during the Middle Miocene to Pliocene (Zhou et al., 2017), even though the plateau and the foreland regions were both deformed during the Paleogene Zapata et al., 2020).
Previous models of the Arizaro region of the Puna Plateau, 150 km north of our study location, also favor a lithospheric dripping scenario, but one in which flexural subsidence formed a large contiguous hinterland basin due to an inferred strong, mafic lower crust (DeCelles, Zandt, et al., 2015;Wang et al., 2015). This crustal strength contrast may be explained by lithologic heterogeneity or inherited structures. It is currently thought that the Arizaro and Antofalla regions are underlain by a contiguous basement terrane, the Antofalla-Arequipa block, onto which multiple episodes of arc magmatism were superimposed during the Phanerozoic (Escayola et al., 2011;Niemeyer et al., 2018;Ramos, 2008;Zimmermann et al., 2014). Although regional geologic studies do not suggest a suture zone between the two regions, we cannot rule out more complex suture geometries or mafic underplating in the lower crust that foundered during the Mio-Pliocene. Arc magmatism can vary in volume and composition along strike (e.g., Cecil et al., 2018;Wei et al., 2021), which could contribute to crustal heterogeneity between the two regions, but more research in the Puna Plateau is necessary to determine this. The two regions also differ in their inherited tectonic structures. A series of continental-scale, NW-SE-trending lineaments cut transversely through the Antofalla region (Richards et al., 2006), but the Arizaro basin is situated between two such lineaments (DeCelles, . Back-arc volcanism (Richards & Villeneuve, 2002;Richards et al., 2006) and crustal shear-wave minima (Bianchi et al., 2013;Ward et al., 2013Ward et al., , 2017 both appear to be correlated with the transverse lineaments. Although the origin of the lineaments is not well constrained, they are thought to be long-lived basement structures that provided conduits for melting and metasomatism throughout the Cenozoic (Richards & Villeneuve, 2002;Richards et al., 2006). This melt and fluid infiltration may have preferentially weakened the lower crust in our study area. Low crustal shear wave velocities also reflect a crust that is weakened by heating or partial melting beneath the Antofalla region (Ward et al., 2017). Lithospheric foundering highlights the geodynamic diversity of major orogens, suggesting an interaction between magmatism, inherited structures, and crustal strength on >10 Myr timescales. More research is required to determine the factors that contributed to such divergent responses to two neighboring lithospheric foundering events.

Conclusions
The role of lithospheric dripping in driving significant crustal deformation has been underappreciated, partly because existing models tend to incorporate strong crusts that resist deformation. We show that lithospheric dripping under a weak crust can cause substantial crustal deformation, including upper-crustal shortening, upper-crustal extension, and lower crustal flow. In the Central Andean Puna Plateau, existing data support the inference that Early-Late Miocene uplift and sedimentation across the SPP are a result of crustal deformation associated with a lithospheric foundering event. Foundering-related deformation persists for several Myr after detachment due to crustal buoyancy contrasts, which can explain the extensional deformation affecting the plateau during the Pliocene to Quaternary. Orogens where lithospheric foundering has been hypothesized likely involved episodes of heterogeneous deformation, uplift, and sedimentation similar to the Puna Plateau.

Data Availability Statement
A detailed description of numerical modeling methods and extended figures and tables are available in Supporting Information S1. Apatite fission track data are available from Carrapa et al. (2005), Stalder et al. (2020), andMcMillan et al. (2022). The Southern Puna earthquake catalog is available from Mulcahy et al. (2014). This work uses deal.II (version 9.2.0) and ASPECT (version 2.4.0-pre, revision 08b6a15), which is open source and can be downloaded at https://aspect.geodynamics.org/. The models of this paper are fully reproducible using the code archived on Zenodo at https://doi.org/10.5281/zenodo.8025049. This work was funded by an NSERC Discovery Grant and Accelerator Supplement to L. Schoenbohm. We thank the Computational Infrastructure for Geodynamics (geodynamics.org) which is funded by the National Science Foundation under award EAR-0949446 and EAR-1550901 for supporting the development of ASPECT. Computations were performed on the Niagara cluster (SciNet/ Digital Research Alliance of Canada).