Controls of the Foreland Deformation Pattern in the Orogen‐Foreland Shortening System: Constraints From High‐Resolution Geodynamic Models

Controls on the deformation pattern (shortening mode and tectonic style) of orogenic forelands during lithospheric shortening remain poorly understood. Here, we use high‐resolution 2D thermomechanical models to demonstrate that orogenic crustal thickness and foreland lithospheric thickness significantly control the shortening mode in the foreland. Pure‐shear shortening occurs when the orogenic crust is not thicker than the foreland crust or thick, but the foreland lithosphere is thin (<70–80 km, as in the Puna foreland case). Conversely, simple‐shear shortening, characterized by foreland underthrusting beneath the orogen, arises when the orogenic crust is much thicker. This thickened crust results in high gravitational potential energy in the orogen, which triggers the migration of deformation to the foreland under further shortening. Our models present fully thick‐skinned, fully thin‐skinned, and intermediate tectonic styles in the foreland. The first tectonics forms in a pure‐shear shortening mode whereas the others require a simple‐shear mode and the presence of thick (>∼4 km) sediments that are mechanically weak (friction coefficient <∼0.05) or weakened rapidly during deformation. The formation of fully thin‐skinned tectonics in thick and weak foreland sediments, as in the Subandean Ranges, requires the strength of the orogenic upper lithosphere to be less than one‐third as strong as that of the foreland upper lithosphere. Our models successfully reproduce foreland deformation patterns in the Central and Southern Andes and the Laramide province.

an orogen may provide first-order control on its basement deformation style. In sediment-starved orogens, such as the Southern Urals case, thick-skinned deformation is mainly located in the orogenic core. When the orogen is sediment-loaded, such as the Swiss Alps, this basement-involved structure appears in the axial zone and its foreland. Moreover, a sudden drop in the mechanical strength of foreland sediments can lead to a shift of the shortening mode from pure-to simple-shear and the formation of thin-skinned structures in the foreland .
However, these studies mainly focused on structural styles of the orogen; the foreland deformation pattern has received less attention. In particular, the exact nature of variations in lithospheric strength and sediment weakening affecting the evolution of foreland deformation is still not well understood. In addition, the question of whether controlled factors from these studies can be applied to explain the deformation patterns in other forelands remains open. The above-cited models also did not explore more details of the foreland deformation features (e.g., the fault direction) due to the lack of necessary numerical resolution. Recent progress in numerical modeling techniques extends this research to higher-resolution lithospheric models, which is the subject of the current study.
The long-term evolution of continental lithospheric strength is primarily controlled by its composition and temperature, which strongly depend on depth, that is, the lithospheric thickness and the crustal thickness (e.g., Cloetingh & Burov, 1996;Ellis, 1988;Kusznir & Park, 1986). The lithosphere can be stronger due to lithospheric thickening and/or crustal thinning ( Figure A1). The composition, fluid content (degree of hydration), magmatism, and thermal/structural inheritance also influence the lithospheric strength (e.g., Burov, 2011;Burov & Watt, 2006;Erdős et al., 2015;Kohlstedt et al., 1995;Mouthereau et al., 2013). For example, the subduction process can weaken the foreland lithosphere in a subduction-dominated orogeny by a high degree of hydration or a hot thermal structure. In this study, we addressed the key (although certainly not all) controlling factors: the thickness of the thermal lithosphere and thickness of the crust. Together, these two factors also automatically determine the partitioning of the lithosphere into the crust and mantle lithosphere, thus also taking into account the effect of the composition and at least partially.
Weakening of foreland sediments can facilitate the initiation of foreland crustal underthrusting below the orogen, thereby promoting the formation of simple-shear shortening . Therefore, the sedimentary strength should also be considered in the development of the foreland deformation pattern. Although the strength of the lithosphere already includes the top sedimentary strength, the latter has a limited effect on the former. On the one hand, although the thickness of the sedimentary cover can reach approximately 8 km or more (Laske et al., 2013), it is still less than one-tenth of a typical continental lithospheric thickness (∼100-200 km thick). Accordingly, changes in sedimentary strength due to thickness variations have little effect on the strength of the entire lithosphere. On the other hand, the shallow frictional brittle strength ( in Equation A6) in the first few kilometers of the crust (1-14 km), which depends highly on pressure rather than compositions (Byerlee, 1978), has less influence on the lithospheric strength than the deep ductile strength ( Figure A1). Thus, changes in the sedimentary strength due to different compositions hardly affect the brittle strength and cause fewer changes in the entire lithospheric strength. In other words, the change in the strength of foreland sediments does not significantly change the strength of the lithosphere. However, it is important for the deformation evolution of the shallow crust of the foreland. Therefore, these two factors should be considered separately regarding the influence on the deformation pattern of the foreland crust.
The friction coefficient of the sediment is another factor other than the thickness that controls its strength (see Appendix A). It has a wide range of values from >0.8 to <0.05, depending on temperature, composition, porefluid pressure, and asperities along the fault surface (Hassani et al., 1997). If sedimentary rocks contain sufficient clay minerals such as montmorillonite or vermiculite (Byerlee, 1978), the friction value can be as low as 0.1. The value can be further decreased to 0.015, which is predicted for subduction channels in some geodynamic models . A reduction in the friction coefficient can decrease the yield strength of the rock, thereby accelerating its failure. The physical nature of potential frictional weakening in foreland sediments remains controversial. This may result from high pore-fluid pressure (lowering the effective confining stress) due to rapid hydrocarbon generation (Cobbold et al. [2004] and reference therein), an increase in precipitation (e.g., Strecker et al., 2007), or compaction under strong compression in the foreland (e.g., . Since we are concerned with crustal-scale deformation in the foreland, the exact origin of the sedimentary friction drop is not discussed here. In this study, we first examine how different factors influence the lithospheric strength of the orogen and its foreland (factors: lithospheric thickness and crustal thickness). We also examine the mechanical strength of foreland sediments (factors: effective friction coefficient of sediments and sedimentary thickness). Then, we systematically investigate how these parameters control the foreland deformation pattern during orogen-foreland shortening. Finally, we apply model results to natural orogen-foreland systems such as the Central-Southern Andes and the Laramide province.

Method and Model Geometry
We use the highly scalable parallel code Lithosphere and Mantle Evolution Model (LaMEM; Kaus et al. [2016]) that solves three geodynamic conservation equations (see Appendix A) to govern material deformation. The model contains two structural domains -the orogen and its foreland -400 km wide and 400 km deep. As we are interested in the deformation of the foreland crust, we plot our modeling results in the zoomed-in area in the top 60 km of the model (dashed gray rectangle in Figure 1) with a horizontal distance between 50 and 330 km. By doing so, we suppose the effect of side boundary conditions on the modeling results in this area minimized (see Figure S1 in Supporting Information S1, showing that the boundary effects on our zoom-in models can be negligible).
The lithospheric thicknesses of the orogen and its foreland in the model vary from 60 to 200 km. Figure 1 shows a 60-km-thick lithospheric strength profile. This is an example of a thin and weak orogenic lithosphere due to lithosphere delamination (e.g., Kay & Kay, 1993). The felsic crust in the foreland is 24 km thick and has a sedimentary cover of varying thickness (0-8 km). Below it, there is a 12-km-thick mafic crust. In contrast, the thickness of the orogenic crust varies between 36 and 70 km. A thick orogenic crust could be produced by tectonic shortening during orogenesis in natural orogens such as the Tibetan Plateau and the Central Andes (e.g., Holt & Figure 1. Initial model geometry with thermal-mechanical boundary conditions. The prescribed compression velocity (V R ) from the right-hand side boundary is balanced by the uniform outflux velocity (V L ) along the left-hand side boundary under the orogenic lithosphere. Orange line is the initial thermal field. The temperature of the lithosphere-asthenosphere boundary (T LAB ) varies between 1324°C and 1380°C, depending on the lithospheric thickness. The crustal thickness in the orogen (H_ oc) varies from 36 to 70 km. The lithospheric thicknesses of the orogen (H_ol) and the foreland (H_fl) vary from 60 to 200 km. The thickness of the foreland sediment (H_se) varies from 0 to 8 km, and the value of its friction coefficient (µ_se) is between 0.5 and 0.02. The white dashed line is the boundary between the orogen and its foreland. Qtz wet , MD dry , and Ol dry in the example of the 60-km-thick lithospheric strength profile indicate wet quartzite, dry Maryland diabase, and dry olivine, respectively. Wallace, 1990;Ramos et al., 2004). The entire model domain has a uniform 500-m-high grid resolution, ensuring that the deformation in such a thin sedimentary layer is tracked correctly.

Material Properties and Boundary Conditions
The material thermomechanical properties are given in Table 1. All materials contain a viscoelastoplastic rheology, where diffusion and dislocation viscous creep regimes are used to mimic ductile deformation. The laboratory-derived flow laws of wet quartzite (Qtz wet ; Gleason & Tullis, 1995), dry Maryland diabase (MD dry ; Mackwell et al. [1998]), and wet/dry olivine (Ol wet /Ol dry ; Hirth & Kohlstedt, 2003) are used for the felsic crust and its sedimentary cover, the mafic crust, and the lithospheric mantle/asthenosphere, respectively. The felsic crust undergoes frictional-plastic strain-softening through a friction coefficient decrease from 0.58 to 0.1 over the accumulated strain of 0.5-1.5, including the friction angle from 30° to 6° and the cohesion from 20 to 1 MPa (Table 1). This follows values used in previous geodynamic models (e.g., Erdős et al., 2015;Sobolev et al., 2006).
The values of thermal parameters are within the range expected for crustal and mantle materials (e.g., Barrionuevo et al., 2021;Sobolev et al., 2006). Radiogenic heat production is 1.0 μW m −3 in the felsic crust and 0.3 μW m −3 in the mafic crust. The thermal conductivity increases from 2.5 W m −1 K −1 in the crust to 3.3 W m −1 K −1 in the mantle. Material density is temperature-dependent (Table 1). The continental felsic crust has a reference density of 2,800 kg m −3 at room temperature to reflect a more felsic (silica-rich) composition than the mafic crust below. The density of the sediments is 300 kg m −3 less than the density of continental felsic rocks at the same temperature. The reference density of the mantle (3,300 kg m −3 ) is consistent with the density of the fertile lithospheric mantle (Poudjom Djomani et al., 2001). Figure 1 shows the initial thermal-mechanical boundary condition. The initial temperature setting of the model is divided into two steps. It first increases linearly with depth from 0°C at the surface to 1328-1380°C at the lithosphere base depending on the lithospheric thickness. It then rises adiabatically to 1460°C at the bottom boundary. The thermal gradient at side boundaries is taken to be zero, which means no horizontal heat flux. We used the "sticky air" top boundary with the free surface stabilization approach (Kaus et al., 2010). This 10-km-thick layer is characterized by low viscosity (10 19 Pa s) and low density (1 kg m −3 ). This boundary condition allows a relatively large integration time step and simulates surface faulting. The boundary condition at the bottom is free-slip. Material flows at a 2 cm/year rate from the right-hand (east) side boundary and out at the left-hand side boundary beneath the orogenic lithosphere to maintain mass balance. The amount of shortening in our models (100 km) is appropriate and reasonable for shortening of the Central Andes over the last 10 Myr (Horton, 2018 Power law exponent a , nf/nl -/4 -/4.7 1/3.5 1/3.5 a Initial temperature-dependent density: ρ P,T = ρ 0 [1-α (T-T 0 )], where ρ 0 is the reference density at temperature T 0 . a Strain-softening in the felsic crust via a decrease in φ and C 0 over the accumulated strain of 0.5-1.5. Sediment is assumed to be initially weak. a Viscous creep includes diffusion (Bf, Ef, Vf, nf) and dislocation (Bl, El, Vl, nl).  , 2006). A higher but reasonable amount of shortening does not change the main results much (see Figure  S1 in Supporting Information S1).

Reference Model
In reference to Model M1, the orogen has the same lithospheric structure as the foreland except for the presence of a 4-km-thick sedimentary layer above the foreland. After 100 km shortening, the felsic crust undergoes pureshear shortening, resulting in distributed crustal thickening and surface uplift ( Figure 2b). Figure 2c shows that the strain rate norm (square root of the second invariant of the deviatoric strain rate) is homogeneously distributed from the surface to the basement at ∼17 km depth. Thus, a fully thick-skinned tectonic style is formed in the foreland. We conducted a series of modeling experiments to systematically investigate how the foreland deformation pattern is affected by changes in the lithospheric structure, crustal structure, and foreland sedimentary strength (Table 2; Figure S2 in Supporting Information S1). Below, we examine the effects of each of the following factors on the deformation style: (a) thickness of the orogenic lithosphere (H_ol), (b) thickness of the orogenic crust (H_ oc), (c) thickness of the foreland lithosphere (H_fl), (d) friction coefficient of foreland sediments (µ_se), (e) thickness of foreland sediments (H_se), and (f) their combinations.

Orogenic Lithospheric Thickness and Orogenic Crustal Thickness
First, we intended to investigate the effect of orogenic lithospheric thickness on the foreland deformation pattern. Geological and geophysical observations indicate that the lithosphere under some active orogens (e.g., the Central Andes) is thin or absent (e.g., Beck & Zandt, 2002;Yuan et al., 2002). This is because the lithospheric mantle, being gravitationally unstable, is susceptible to partial removal via Rayleigh-Taylor-type instability (Molnar & Houseman, 2004) or complete removal by delamination (e.g., Bird, 1979;Le Pourhiet et al., 2006). In Model M2 (Figure 3a), the orogen has a 60-km-thick lithosphere and is weaker than the foreland in which the lithosphere is 100 km thick. The model result shows that the compressional deformation is localized within the orogen after 100 km pure-shear shortening. Simultaneously, a fully thick-skinned structure is formed in the foreland. In contrast, when the initial lithosphere of the orogen is thicker and therefore stronger than that of the foreland (as in Model M3 in Table 2), the foreland deformation pattern remains unchanged. Therefore, in the models where only the orogenic lithospheric thickness changes while the crustal thicknesses in the orogen and its foreland remain the same, shortening of the foreland crust is in pure-shear mode accompanied by a fully thick-skinned tectonic style.
When the initial crust of the orogen thickens to 60 km, the foreland crust underthrusts beneath the orogen regardless of the thickness of the orogenic lithosphere within the range of parameters considered (Table 2), which is interpreted as a simple-shear shortening mode. In this mode, if the contribution of the thin-skinned deformation to the total foreland crustal deformation is less than one-tenth, then we consider this tectonic style to be thick-skinned dominated (e.g., M4-M6). The amount of simple-shear deformation appears to be greater in Model M4 (Figure 3b). It has a thinner orogenic lithosphere than in Model M6 (Figure 3c). In both models, a pronounced deep shear zone is produced between the upper and lower crust in the foreland.

Foreland Lithospheric Thickness
Here, we tested the effect of the foreland lithospheric strength on the deformation style by changing the foreland lithospheric thickness. In contrast, the initial crustal thicknesses in the foreland and the orogen were fixed. Unlike in the mountain belts, the foreland lithosphere in the craton can be thicker than 150 km. For example, the thermal lithosphere is > 180 km thick under some foreland regions of the southwestern Canadian craton (Currie, 2016). In Model M8 with a 200-km-thick foreland cratonic lithosphere, most of the shortening is concentrated in the orogenic crust, resulting in crustal buckling and surface uplift (Figure 3d). A fully thick-skinned structure is formed near the orogen-foreland boundary. As expected, the amount of foreland deformation decreases with the thickening of the foreland lithosphere (e.g., comparing M1 and M8). Note. H_ol: thickness of the orogenic lithosphere, H_fl: thickness of the foreland lithosphere, H_oc: thickness of the orogenic crust, H_se: thickness of foreland sediments, µ_se: friction coefficient of foreland sediments; S. mode: shortening mode, S-1: pure-shear, S-2: simple-shear; T. style: tectonic style, T-1: fully thick-skinned, T-2: thick-skinned dominated, T-3: thin-and thick-skinned mixed, T-4: fully thin-skinned.

Table 2
The List of the Orogen-Foreland Shortening Models

Foreland Sedimentary Strength
The foreland sedimentary strength (coefficient of friction and its thickness) is also important for the foreland deformation pattern. Here, we varied the friction coefficient values of foreland sediments in the range of 0.58-0.02 (M1, M9-M10). This range is consistent with that of previous geodynamic models (e.g., Sobolev et al., 2006). The foreland deformation in Model M9 is no longer homogenous as in the reference model; pronounced thrust faults are produced in the middle part of the foreland (Figure 3e). When the friction coefficient of sediment is further reduced, and its thickness increases (M10), the magnitude of deformation in the foreland increases, and the fault system become more complicated. However, the deformation pattern retains pure-shear shortening and fully thick-skinned tectonics. Regarding the factor of the thickness of foreland sediments, our model results show that the deformation pattern does not change if only the sedimentary thickness is increased, but less and deeper faulting occurs in the foreland crust ( Figure S2c in Supporting Information S1).  Table 2 after 100 km shortening, showing a-e) effects of individual factors and f-j) effects of multiple factors. Foreland sediments (the red part in the structure bar of the foreland lithosphere) are considered initially weak when their thickness is greater than 4 km and their friction coefficient is not higher than 0.05. The black solid line with two arrows represents the tectonic style of thin-skinned in the foreland.

Effects of Multiple Factors
None of the above models shows a wide thin-skinned thrust zone in the foreland. Here, we present models combining multiple factors considered above (Figures 3f-3j). All of these models have a 4-km-thick sedimentary layer in the foreland with a friction coefficient of 0.05 (we term "weak foreland sedimentary layer", i.e., the red area in the lithospheric structure bar in Figure 3). As we will see later, weak foreland sediments result in two additional tectonic styles thin-and thick-skinned mixed and fully thin-skinned. We deem the tectonic style to be mixed if its thin-skinned thrust zone is significantly wider than the zone in thick-skinned dominated tectonics (e.g., Figure 3i).
Weakening of foreland sediments in most models facilitates the underthrusting of the foreland beneath the orogen and promotes the formation of thin-and thick-skinned mixed or fully thin-skinned tectonics. The formation of the latter further requires a relatively thicker crust and/or thin lithosphere in the orogen (M12, M14, M17, M20). Foreland weak sediments can also switch the shortening mode from pure-shear to simple-shear (Figures 3a  and 3f) when the crust in the orogen is thin but its lithosphere is thinner than the foreland lithosphere. This switch does not occur if the orogenic lithosphere is thicker (e.g., compare M3 and M7 with M13 and M15) or if the thicker foreland is in the cratonic area (Figures 3d and 3h). Additionally, these combined models show that large foreland underthrusting and mid-crustal viscous flow lead to crustal thickening and surface uplift in the orogen (Figures 3g and 3j).

Lithospheric Strength Analysis
We calculated the initial integrated lithospheric strength of the orogen and its foreland and the strength ratio between them. The integrated strength is estimated through the integration of the yield strength envelope (e.g., Burov, 2011;Tesauro et al., 2013). Since the strength of the relatively thin sedimentary layer has little effect on the lithospheric strength, we neglected the strength change caused by the weakening of foreland sediments during the calculation. More details about the calculation are presented in the Appendix A.
As we show below, foreland deformation styles are first-order controlled by the difference in lithospheric strength between the orogen and the foreland (Figure 4). We note, however, that the difference in the integrated strength of the entire lithosphere between the orogen and the foreland does not explain all model results. For example, the entire lithospheric strength of the orogen in Model M18, including a 150-km-thick lithosphere and a 60-km-thick crust, is higher than that in Model M21 with an 80-km-thick lithosphere and a 36-km-thick crust in the orogen ( Figure A1, Figures S1 in Supporting Information S1). Model M18 behaves in simple-shear shortening with thin-and thick-skinned mixed structures in the foreland ( Figure 3i). As expected, if the model has a thinner and thus weaker orogenic lithosphere, it shows further underthrusting of the foreland beneath the orogen and a larger amount of thin-skinned deformation in the foreland (e.g., compare M14 with M18). However, the model behavior of M21 contradicts this view, where the tectonic type is thick-skinned dominated with a narrow thinskinned wedge zone on the edge of the foreland ( Figure S2e in Supporting Information S1).
Therefore, we considered the strength difference of the upper part of the lithosphere between the orogen and its foreland, which controls the foreland deformation pattern better than the strength difference of the entire lithosphere ( Figure 4). With this new definition of the upper lithospheric strength, Model M21 has a stronger upper orogenic lithosphere than Model M18, and therefore its strength ratio is higher. As a result, M21 shows less thinskinned deformation in the foreland.
If the upper lithospheric strength in the orogen and its foreland are similar (strength ratio ∼0.8-1.3 in Figure 4), the foreland (and the orogen) should deform in a pure-shear shortening mode accompanied by thick-skinned deformation. Less obvious is the foreland simple-shear mode and thin-skinned tectonics at a low strength ratio, that is, when the orogenic lithosphere is much weaker than the foreland lithosphere. In this case, the intuitive scenario is the localization of shortening in the weak orogen rather than the foreland. However, the strong foreland in our models behaves in different deformation patterns. We infer that in addition to the lithospheric strength mentioned above, the gravitational potential energy (GPE) of the orogen also contributes to the foreland deformation pattern.
Generally, the compressive force driving orogenic shortening (i.e., mountain building) causes the thickening of the orogenic crust. During shortening, the force works against two main resistive forces: mechanical strength (discussed in this study) and gravity (e.g., Molnar & Lyon-Caen, 1988). The work against gravity creates gravitational potential energy. The GPE per unit surface of the Earth area in the orogen increases with crustal thickening. Thus, shortening the orogen further requires an increasingly larger amount of work from the driving force to overcome the increasing GPE. When the force can no longer supply the energy needed to elevate the orogen higher, the mountain range is likely to grow laterally in width instead of increasing in height and crustal thickness (Molnar & Lyon-Caen, 1988). Consequently, when the orogen grows laterally, the work done by the specified driving force will be used to deform the orogenic edge and its foreland, even if the orogenic lithosphere is much weaker than the foreland lithosphere. In this scenario, the foreland lithosphere can underthrust beneath the edge of the orogen, that is, the foreland shortening mode is simple-shear ( Figure 4). If there is a thick layer of mechanically weak sediments in the foreland, shear deformation is localized in the sedimentary layer and the foreland tectonic style is thin-skinned (Figures 4b and 4d). In this study, we treat the role of GPE as a reasonable qualitative assumption without testing its effect on lithospheric strength because the GPE of the orogen is in tur n controlled by its crustal thickness and lithospheric thickness.

Structural Controls on the Shortening Mode and Tectonic Style in the Foreland
Our model results demonstrate that the variation in orogenic strength caused by the change in orogenic crustal thickness has a critical effect on controlling the shortening mode. The pure-shear mode develops in the models with little difference in the crustal thickness between the orogen and the foreland, while the thickened orogenic crust is required to switch from pure-shear to simple-shear (Figure 4). The thickened orogenic crust causes the initially high GPE of the orogen and low strength of the orogenic upper lithosphere. This high GPE forces the shortening shift to the foreland. The thick and weak orogenic crust allows the strong foreland lithosphere to intrude into it in simple-shear mode easily. Our models show that the other four individual factors (H_ol, H_fl, µ_sed and H_sed) have little effect on the transition of the shortening mode with one exception. That is the case (the dashed rectangle in Figure S2b in Supporting Information S1) when the orogenic crust is much thicker (high GPE) than the foreland crust and the foreland lithosphere is thin, showing a pure-shear shortening mode in the foreland.
Our models show that the significantly lower strength of the upper lithosphere in the orogen than in the foreland (strength ratio < ∼0.7) and the presence of thick and weak foreland sediments are responsible for the thin-skinned tectonics in the foreland. The absence of these conditions results in the tectonic style of fully thick-skinned or thick-skinned dominated. Furthermore, the condition of thick and weak foreland sediments generally intensifies simple-shear shortening by making underthrusting easier and thus broadening the thin-skinned thrust zone. When the orogenic crust is thick and the foreland lithosphere is thin, this condition can switch the shortening mode in the foreland from pure-shear to simple-shear.

Applications to Natural Orogen-Foreland Shortening Systems
Here, we compare our model inferences to the Central and Southern Andes and the Laramide orogeny and provide a first-order fit of the foreland deformation pattern to these natural shortening systems. We will look more specifically at the Altiplano-Puna plateau-foreland profile (Figures 5b and 5c), the Frontal Cordillera-Precordillera-Sierras Pampeanas profile (Figure 5d), and a more conceptual cross-section through the Colorado Plateau and Southern Rocky Mountain foreland (Figures 5e and 5f).

Altiplano-Puna Plateau
In the Central Andes, the Altiplano-Puna Plateau was formed with N-S oriented deformation diversity, including a broad wedge-shaped thin-skinned thrust belt in the Interandean-Subandean zone and a thick-skinned structure in the Santa Barbara System (Figure 5a). The lithosphere under the plateau is very thin, but the upper felsic crust is as thick as 50-70 km (e.g., Ibarra et al., 2019;Tassara et al., 2006). This inherited thin lithosphere is suggested to result from lithospheric delamination, which occurred during Cenozoic shortening (e.g., Beck & Zandt, 2002;Kay & Coira, 2009;Kay & Kay, 1993;. The Puna Plateau and its foreland area have a higher seismic attenuation, implying a hotter and thinner lithosphere than the northern Altiplano part (Whitman et al., 1996). Paleozoic and Mesozoic sediments were abundantly deposited in the Subandean zone but pinched out southward to the Santa Barbara system (e.g., Allmendinger & Gubbels, 1996;Pearson et al., 2013). The local wet conditions in the foreland since the late Cenozoic (Strecker et al., 2007) indicate that abundant fluids are stored in these ancient sediments and may weaken them by increasing their pore fluid pressure.
We applied these observations to the model of the Central Andes. In the models (Figures 5b and 5c), the thickness of the orogenic crust under the Altiplano-Puna Plateau is 60 km. An additional 10-km-thick lithospheric mantle is attached to the Altiplano crust. The orogenic lithosphere under the Puna Plateau only contains a thick crust due to mantle lithospheric delamination. The lithosphere of the Puna foreland in the model is 70 km thick and 30 km thinner than the Altiplano foreland lithosphere. In agreement with observations, the weak sedimentary layer in the model covers only the north Altiplano foreland crust (Figure 5b). Model results clearly show that the simpleshear mode with a fully thin-skinned thrust belt and the pure-shear mode with a fully thick-skinned structure are formed in the Altiplano foreland and the Puna foreland, respectively. Our models support and specify the results of previous relatively low-resolution modeling studies (e.g.,  and reproduce observed east-dipping reverse faults in the foreland edge in both cases.

Precordillera-Sierras Pampeanas Region
The Sierras Pampeanas province, located on the eastern side of the Precordillera thin-skinned thrust belts, is a modern analog of the thick-skinned deformation of the Laramide province (Jordan & Allmendinger, 1986). The tectonic style of the Precordillera-Sierras Pampeanas foreland region, adjacent to the Frontal Cordillera, can be broadly considered a thin-and thick-skinned mixed structure (Figure 5a). The oceanic flat slab below the Frontal Cordillera stays at 100 km depth, and thus, the orogenic lithosphere of the Frontal Cordillera may be less than 100 km thick (e.g., Jordan et al., 1983;Ramos & Folguera, 2009). The lithospheric thickness increases eastward and is more than 20 km thicker in the Sierras Pampeanas foreland. Crustal thickness exceeds 60 km beneath the Frontal Cordillera and rapidly decreases eastward to less than 40 km below its foreland (e.g., Perarnau et al., 2012;  Lacombe and Bellahsen (2016) and Yonkee and Weil (2015), and its modeled foreland deformation pattern. Ramos et al., 2004). Furthermore, there are abundant Paleozoic sedimentary rocks in the Precordillera, whereas only a small amount of Cenozoic sediments covers the Sierras Pampeanas (e.g., Meeβen et al., 2018;Ramos et al., 2004).
Unlike the 30°-dipping subducted slab in the Central Andean case, the southern Argentine Andean case slab is nearly horizontal (Gutscher et al., 2018;Jordan et al., 1983). Slab flattening can enhance the stress transmission from the subducting plate into the overlying plate by increasing the degree of plate coupling (e.g., Lacombe & Bellahsen, 2016), thus promoting plateau-foreland shortening, which may contribute to the development of thick-skinned tectonics. Note, however, that in the cases of the Sierras Pampeanas and the Laramide below, we do not introduce the factor of flat-slab subduction; therefore, our models do not exactly reproduce the amount of shortening and high topography of these two provinces.
The model constrained by these observations includes a thin and weak orogenic lithosphere 30 km thinner than the foreland lithosphere. Crustal thickness is greater than 60 km in the orogen and decreases to ∼40 km in the foreland. The model result (Figure 5d) indicates that simple-shear shortening occurs in the foreland, accompanied by mixed tectonics consisting of thin-skinned thrust at the foreland edge (Precordillera) and thick-skinned structure behind (Sierras Pampeanas). Note that the weak sedimentary layer is located through the entire foreland area in our model. Although it has little influence on the tectonic style of the Sierras Pampeanas, it is necessary to consider the difference in sedimentary thickness between the Precordillera and Sierras Pampeanas in future studies.

Laramide Province
The Laramide province (i.e., the Rocky Mountain foreland adjacent to the Colorado Plateau) is a widely thickskinned deformation zone that developed more than 1,000 km inboard the plate margin (e.g., Bird, 1984;Erslev, 2013;Lacombe & Bellahsen, 2016;Saleeby, 2003;Yonkee & Weil, 2015). This province sustained more than 100 km pure-shear shortening, which contrasts strongly with minor deformation of the Colorado Plateau (e.g., Bird, 1984;Flowers et al., 2008;Humphreys, 2009;Spencer, 1996). Dynamic processes that propagate deformation across the strong and broad plateau far into the foreland and produce thick-skinned tectonics in the case of the Laramide orogeny are still largely debated.
One fashionable possibility is that the formation of the Laramide province is suggested to be the result of slab flattening of the Farallon plate. In particular, this process enhanced interplate coupling along the base of the cratonic lithosphere root, hence efficient stress transmission from the Farallon plate into the North American plate to deform the plateau-foreland system (e.g., Axen et al., 2018;Bird, 1984). Furthermore, flat-slab subduction likely changed the strength of the continental lithospheric mantle. For instance, a cold slab can cool the above basal lithospheric mantle, which favors increased strength and stress transfer far into the foreland. In contrast, the lithospheric mantle can also be weakened due to the effects of basal lithospheric mantle removal by flat-slab subduction (e.g., Axen et al., 2018;Bird, 1984;Liu & Currie, 2016), hydration from dewatering of the underlying flat slab and heating by a magmatic ascent (Humphreys et al., 2003), and/or thermal inheritance from the preorogenic extension (Marshak et al., 2000). Lithospheric mantle weakening may allow shortening to occur in the deep mantle beneath the southern Rocky Mountains. Together with enhanced stress transfer, this process possibly promotes crustal shortening and leads to thick-skinned deformation within the foreland.
In addition to flat slab subduction, crustal/lithospheric buckling has been considered another possible mechanism for propagating and accommodating deformation in the Laramide foreland (e.g., Erslev [1993]; Lacombe & Bellahsen [2016]; Tikoff & Maxson [2001] and reference therein). For instance, Lacombe and Bellahsen (2016) emphasize that thick-skinned tectonics in the orogenic foreland is favored by the occurrence of a ductile middle or lower crust of a young, and hot lithosphere, hence enabling crust-mantle decoupling. Depending on its composition -felsic or mafic granulites -the middle or lower crust may have been either moderately weak with the potential concentration of ductile flow along deep décollements or strong with potential for lithospheric buckling (Yonkee & Weil, 2015). Overall, intervening in specific boundary conditions such as flat slab subduction, together with structural crustal inheritance and possible mantle weakening, may provide a sophisticated explanation for intraplate basement-involved shortening in the Laramide setting.
As the deformation did not regularly propagate in a classical 'in sequence', foreland-ward way from the former Sevier orogen to the Laramide orogen, individual basement-cored deformation zones in the Laramide province may have developed spatially and temporally in a rather complex sequence (e.g., Crowley et al., 2002;Lacombe & Bellahsen, 2016). Since we are concentrated on the foreland deformation pattern during the Laramide orogeny, we simply developed a conceptual plateau-foreland shortening model constrained by observations of SW-NE tectonic transect through the Colorado Plateau and Southern Rocky Mountain foreland. This does not include the Sevier orogeny (Figure 5e). Although both western Farallon flat slab subduction and eastern intraplate shortening between the Colorado Plateau and the Rocky Mountain foreland can occur during Laramide deformation, we focus on the latter event. The subduction process is not integrated into the Laramide shortening model. Alternatively, we suppose that the presumptive flat-slab subduction on the left boundary prevents the leftward motion of the plateau. Hence, we close the left boundary above the orogenic lithosphere, resulting in a high degree of coupling between the plateau and its foreland.
In this transect, the Colorado Plateau and nearby Rocky Mountain foreland presumably involved a cool and thick lithosphere during the Laramide orogeny. Xenolith-based observations estimate the lithospheric thickness of the Colorado Plateau to be more than 150 km due to its underlying cold, refractory mantle root (e.g., Li et al., 2008;Smith & Griffin, 2005). Previous numerical studies of the flat slab subduction suggest that the Colorado Plateau may be thicker and thus stronger than its foreland cratonic lithosphere due to its deep cratonic root (e.g., Liu & Currie, 2016;O'Driscoll et al., 2009). The foreland was formerly part of a continental platform with an approximately 33-km-thick crust before the Laramide orogeny (Bird, 1984). The lower crust is cool, viscous, and largely intact beneath the Colorado Plateau (Humphreys et al., 2003). Lithostratigraphic columns of Laramide sedimentary successions in depocenters of key Laramide basins show that the thickness of the sedimentary cover is ∼1-2 km (Dickinson et al., 1988).
To develop the model, we used a similar setup to the previous geodynamic models of the Laramide (Liu & Currie, 2016) and further constrained it with the above observations. In this model, the plateau lithosphere is 240 km thick, 40 km thicker than the foreland lithosphere. The foreland crust is 10 km thinner than the orogenic crust and includes a 2-km-thick sedimentary layer. When the strength of the upper lithosphere of the orogen is slightly greater than that of the foreland (the hollow star of the Laramide-R.M. case in Figure 4c), the foreland is subjected to pure-shear shortening with fully thick-skinned tectonics (Figure 3f), and there is minor deformation in the plateau. Foreland deformation is mainly accommodated in the felsic upper-middle crust, potentially implying decoupling between the felsic upper-middle crust and mafic lower crust and lithospheric mantle. Our model results support the mechanism of lithospheric buckling in Laramide deformation.
Note that we have not attempted to provide a thorough review of the Andean/Laramide orogeny. Rather, we have attempted to demonstrate that the foreland deformation pattern of the Andean/Laramide orogeny is consistent with simplified orogen-foreland shortening models. The very fine internal structure of the deformed sediments is not visible in our models and is modeled as a zone with a finite strain greater than 1. This is because our models did not employ a deformed mesh used in Erdős et al. (2015) and Jammes and Huismans (2012), although the resolution of our models is sufficient. We have addressed only the contrast in the lithospheric strength between the orogen and foreland and the strength of the foreland sediment within these shortening models. Other parameters (e.g., the rate and amount of shortening, subduction dynamics, and thermal/structural inheritance) have not been addressed here but are necessary to be considered in future comprehensive case studies.

Conclusions
With high-resolution thermomechanical numerical models, we systematically examine the effects of the lithospheric structure and foreland sedimentary strength on the foreland deformation pattern subjected to tectonic shortening.
We find that three factors significantly control the shortening mode (pure-shear or simple-shear) and the tectonic style (thick-skinned or thin-skinned): (a) the strength difference in the upper lithosphere between the orogen and its foreland, rather than the difference in the entire lithospheric strength between them, (b) GPE of the orogen that is in turn controlled by its crustal thickness and lithospheric thickness, and (c) the strength and thickness of the deforming foreland sediments.
If the strength of the orogenic upper lithosphere is higher or similar to that of the foreland upper lithosphere (strength ratio >∼0.8) and the orogenic crust is not much thicker than the foreland crust (relatively low GPE of the orogen), pure-shear shortening develops in the foreland.
If the strength of the orogenic upper lithosphere is significantly lower than that of the foreland upper lithosphere (strength ratio <∼0.7) and the orogenic crust is much thicker than the foreland crust (>50 km causing relatively high GPE of the orogen), the foreland undergoes simple-shear shortening.
Fully thin-skinned or thin-and thick-skinned mixed tectonic styles can develop in the foreland only if thick (>∼4 km) and mechanically weak (friction coefficient < ∼0.05) sediments are present in the simple-shear shortening mode. Furthermore, the most pronounced fully thin-skinned tectonics develops in the thick and weak foreland sedimentary layer when the strength of the orogenic upper lithosphere is much lower than that of the foreland upper lithosphere (strength ratio <0.3-0.4; Altiplano-Subandean ranges case).
Our high-resolution orogen-foreland shortening models successfully reproduce foreland deformation patterns in the Central and Southern Andes in South America during the Neogene and Laramide province in North America during the Late Cretaceous to Paleocene.

Appendix A: Geodynamic Governing Equations and Yield Strength Envelope
Material deformation is governed by solving the coupled system of momentum (Equation A1), mass (Equation A2), and energy (Equation A3) conservation equations below: where i, j represent spatial directions following Einstein summation convention, x i,j are the Cartesian coordinates, ij is the deviatoric stress tensor, P is pressure, ρ is the density, g i is the gravitational acceleration vector, and are components of the velocity, D/Dt is the material time derivative, C p is specific heat, k is thermal conductivity, A is the radiogenic heat production, and ̇ v ij ,̇ p ij are viscous and plastic strain-rate deviators, respectively. Repeated indices imply summation. These basic geodynamic equations are solved assuming plane strain, incompressibility, and neglecting thermal diffusion.
The material behaves the frictional-plastic deformation when the deviatoric stress exceeds the plastic yield stress ( Y ), which follows a pressure-dependent Drucker-Prager yield criterion: where φ is the internal friction angle and C 0 is the cohesion. We assume the friction coefficient = ( ) . Below this yield stress, materials deform viscously with an effective viscosity ( eff ) given by: where h is the lithospheric thickness, 1 and 3 are the maximum and minimum principal stress components, respectively. Figure A1 shows initial strength envelopes of the lithosphere with different structures. There are two different types in the envelope: the frictional brittle strength ( B ; the purple line in Figure A1) and the ductile strength ( D ; dashed colored curves in Figure A1). The brittle strength is estimated by Byerlee's law (Byerlee, 1978) and a function of pressure-independent rock types in a compressional environment: B = ∫ . The initial reference strain rate ( ref ) is 10 −16 s −1 . The viscous parameters are corresponding to the dislocation creep mechanism. Figure A1. The list of initial lithospheric strength curves with different initial lithospheric structures (60-200 km) and crustal structures (36-60 km), showing lithospheric strengths of the orogen and its foreland for each of the aforementioned model in Table 2. For example, M1-M5: Foreland, means that models M1-M5 contain an initial 100-km-thick lithosphere in the foreland. M2, M12: Orogen, means that in M2 and M12, the initial thickness of the orogenic lithosphere is 60 km and its crust is 36 km thick.

Data Availability Statement
The open-source code LaMEM can be found at https://bitbucket.org/bkaus/LaMEM. For this paper, we used the master branch and commit version b58966b. Computations were performed with resources provided by the North-German Supercomputing Alliance. Input files used to produce the model results are available at https:// doi.org/10.5281/zenodo.5963016.