Kinematic reconstruction of the Bolivian orocline

Determining the relationship between crustal shortening and surface uplift of orogenic plateaus remains a frontier problem in continental geodynamics. The second-largest modern orogenic plateau in the central Andes has been of particular interest because there is evidence for protracted deformation and crustal shortening since at least 40 Ma (Elger et al., 2005; Horton, 2005; McQuarrie et al., 2005; Oncken et al., 2006), but also recent surface uplift as young as 10 Ma that would postdate significant shortening (e.g., Gregory-Wodzicki, 2000; Garzione et al., 2008; Hoke and Garzione, 2008). From southern Peru to northern Argentina (∼14°S–22°S), the Andes feature a high shortening retroarc fold-thrust belt (123–330 km, 37%–40%; Gotberg et al., 2010; McQuarrie, 2002) with peak elevations to ∼5–6 km that form the eastern boundary of the ∼3–4-km-high Altiplano plateau (Fig. 1A). In addition, the region is characterized by a prominent bend in both the continental margin and cordillera from Peru to Bolivia that is referred to as the Bolivian orocline (Carey, 1955) (Fig. 1). Work in the Altiplano–Bolivian orocline area has developed two differing bodies of data that argue for (1) protracted isostatic surface uplift due to tectonic shortening and crustal thickening (e.g., Isacks, 1988) or (2) rapid surface uplift in response to removal of the lower lithosphere (e.g., Garzione et al., 2008). In either case, orogenic shortening should be able to account for both the modern crustal thickness as well as additional material that may have been removed. Available shortening estimates vary from fully accounting for modern crustal thickness (Müller et al., 2002; McQuarrie, 2002; Elger et al., 2005; McQuarrie et al., 2005) to predicting crustal deficits (Kley and Monaldi, 1998; Gotberg et al., 2010). However, the application of these two-dimensional (2-D) plane-strain shortening estimates to crustal thickness calculations is complicated by the kinematic system of the Bolivian orocline. Lithospheric thickening and delamination events, similar to scenarios proposed for the Altiplano, have been related to the 3-D kinematic system associated with orocline formation (Gutiérrez-Alonso et al., 2004, 2011).

Unraveling the complex 3-D kinematics of the Bolivian orocline requires accounting for material displacements parallel to the orogenic trend (Kley, 1999; Hindle et al., 2005; Arriagada et al., 2008) and vertical axis rotations (Arriagada et al., 2006; Roperch et al., 2006) that are not resolved in 2-D cross sections. Internal strain is normally a critical part of the 3-D kinematic system in curved orogens (Yonkee and Weil, 2010), but has been shown to be negligible in the Bolivian orocline (Eichelberger and McQuarrie, 2015). Considering these factors in addition to age constraints is critical to accurately determining the influence of orocline deformation on the evolution of the crust and mantle lithosphere at the Altiplano. To achieve this, we present new high-resolution reconstructions of the Bolivian orocline with the goal of developing a more accurate kinematic model for the central Andes and Andean Plateau.

The reconstructions presented here focus on defining the 3-D kinematics of the Bolivian orocline and investigating how out-of-plane displacements have influenced crustal shortening at the central Andes and oroclines in general. Vertical axis rotation and strike-slip faulting have been shown to be the main mechanisms accommodating out-of-plane motion during the development of orogenic curvature (e.g., Kley, 1999; Arriagada et al., 2008). While the location and sense of offset on faults accommodating strike-slip displacements are now fairly well defined (Dewey and Lamb, 1992; Funning et al., 2005; Eichelberger et al., 2013), there are no direct controls on offset magnitude (Eichelberger et al., 2013). Permissible fault-parallel displacements can be estimated from reconstructions, but the reliability of these estimates is largely dependent on how well constrained the surrounding kinematics are. These reconstructions incorporate variable shortening magnitudes and directions; a range of possible limb rotations from global positioning system (GPS) and paleomagnetic data; and timing constraints from geology and thermochronology into a tightly constrained kinematic framework. Based on the results of these models, we can determine the collective effect of rotation and fault parallel motion on shortening at the Bolivian orocline.

Out-of-plane motion has been shown to be a substantial factor in the development of orogenic curvature (e.g., Hindle and Burkhard, 1999; Yonkee and Weil, 2010). In the central Andes, regional structural trends and paleomagnetic data imply that the Bolivian orocline is no exception. A northwest-southeast–trending structural fabric extending from southern Peru to central Bolivia defines the northern limb of the orocline (Fig. 1B). At the orocline axis in central Bolivia, this transitions to a north-south structural fabric that defines the southern orocline limb southward into Argentina and Chile. Based on geologic map patterns, the principle shortening directions of the limbs predict extension at the orocline core, but none is observed (Figs. 1B and 2A) (Kley, 1999; Eichelberger et al., 2013). Instead, paleomagnetic data and mapped strike-slip displacements indicate that the divergence of shortening directions in the limbs is accommodated by out-of-plane displacement (Kley, 1999). Paleomagnetic data from across the central Andes characterize a regional pattern of counterclockwise rotations in the northern limb from Peru to northern Bolivia and clockwise rotations in the southern limb from southern Bolivia into northern Chile and Argentina (Fig. 3) (Arriagada et al., 2006; Barke et al., 2007; Roperch et al., 2006, 2011). Geodetic data indicate that orogenic curvature continues to evolve today with the potential for 4°–10° of limb rotation over the past 10 m.y. (Allmendinger et al., 2005). In addition to vertical axis rotations, earthquake focal mechanisms and geologic mapping indicate fault-parallel displacements on the Cochabamba fault (Dewey and Lamb, 1992) and Rio Novillero fault (Dewey and Lamb, 1992; Eichelberger et al., 2013) at the orocline core (Fig. 1). In prior reconstructions of the Bolivian orocline, fault-parallel displacements are crucial in restoring shortening in the absence of significant extension (Kley, 1999; Eichelberger et al., 2013).

The Carey (1955) classification of the central Andean curvature as an orocline carried the kinematic implication that the orogen was initially straight and underwent subsequent limb rotation and regional bending. Later work suggested that curvature developed due to a shortening gradient, increasing from the orocline limbs toward the axis in central Bolivia (Isacks, 1988). However, the differences in Bolivian retroarc shortening are insufficient to fully account for paleomagnetic and GPS rotations, which suggests that curved slips paths or regional bending form a significant component of rotation (Eichelberger et al., 2013). As such, the modern orogenic curvature cannot be attributed to a shortening gradient alone, but the temporal overlap in vertical axis rotation and shortening indicate curvature developed contemporaneously with fold-thrust belt deformation.

Given the complexity of orocline kinematics, map-view reconstructions have long been recognized as the best tool for resolving shortening normal to structural trends as well as orogen-parallel motion (e.g., Laubscher, 1965). Where data and geologic constraints are sparse or nonexistent, the requirement of kinematic compatibility can provide predictions on displacement directions and magnitudes necessary for a self-consistent reconstruction (McQuarrie and Wernicke, 2005). Early reconstructions leveraged kinematic compatibility to predicted ∼350 km of shortening at the orocline axis for limb rotations of 5°–10° (Kley, 1999). However, the kinematic resolution of the model was limited by sparse geologic mapping and paleomagnetic data as well as shortening estimates that underpredicted crustal thickness (Kley and Monaldi, 1998). Models constrained by paleomagnetic data (rather than shortening estimates) predicted ∼400 km of shortening at the orocline axis to restore forearc rotations that average ∼30° (Arriagada et al., 2008). These models focused on regional orocline bending rather than the specific kinematics of orocline deformation. Using updated structural geometry for the orocline core, simplified reconstructions predicted ∼100 km of orogen-parallel strike-slip displacements with ∼14° limb rotations (Eichelberger et al., 2013) but were limited in spatial extent.

None of the prior reconstructions fully incorporate both shortening estimates and paleomagnetic data at the structural resolution available in geologic maps. Mapping from the region (Kley, 1996; McQuarrie, 2002; Müller et al., 2002; Elger et al., 2005; McQuarrie et al., 2008; Eichelberger et al., 2013) has provided detailed control on variations in displacement direction and shortening magnitude that were unavailable or greatly simplified in previously published reconstructions. In general, these reconstructions also restore deformation in a single time step, but sufficient thermochronologic and geologic data sets now exist (Sempere et al., 1990; DeCelles and Horton, 2003; Horton and DeCelles, 1997; Horton, 2005; Uba et al., 2005, 2006, 2009; Gillis et al., 2006; Barnes et al., 2006, 2008, 2012; McQuarrie et al., 2008; Eichelberger et al., 2013) that constrain variations in deformation timing across the fold-thrust belt. The temporal component of deformation has been a missing aspect of prior reconstructions but is a critical variable in relating crustal thickening driven by orocline deformation to the timing of surface uplift.

To reconstruct orocline kinematics in a geographic context, this study utilizes an ArcGIS (geographic information system) based map-view reconstruction (Fig. 4A). The model covers the fold-thrust belt region of the Bolivian orocline and sequentially restores deformation from ca. 50 Ma to present in 5 m.y. increments. The restoration is driven using a Microsoft VBA (Visual Basic for Applications) script that calculates the positions of polygon centroids, displacing them at each time step based on vectors defined in an associated dBASE stable (after McQuarrie and Wernicke, 2005). The model accounts for vertical axis rotations by prescribing rotations about a polygon centroid, but for regional rotations each polygon’s restoration path must be manually adjusted to reflect the regional change in displacement direction.

The viability of the proposed restoration paths can be assessed based on the self-consistency of the kinematic model through time (McQuarrie and Wernicke, 2005). Polygons that restore to overlapping positions indicate regions that violate strain compatibility due to defined restoration paths that are kinematically incompatible (McQuarrie and Wernicke, 2005). In the Bolivian Andes, deformation is principally characterized at both outcrop (McQuarrie and Davis, 2002) and regional scales (e.g., Sheffels, 1990; Kley, 1996; Lamb and Hoke; 1997; McQuarrie, 2002) by structures accommodating horizontal shortening perpendicular to structural trend. Balanced cross sections account for this plane-strain shortening; thus restoring Andean shortening estimates in map view requires the displacement direction be parallel to the line of section orientation and perpendicular to structural trend. Consequently, the restoration vectors in our reconstructions are initially defined by the magnitude of shortening directed normal to structural trend (Fig. 2A). Given the similarity of shortening estimates in the Bolivian Andes, kinematic incompatibilities largely arise from regional changes in structural orientation that produce nonparallel restoration paths (Fig. 2B). Here we interpret the locations of kinematic incompatibilities to indicate regions where plane-strain deformation constraints are insufficient and additional components of out-of-plane motion are required (overlap in Fig. 2B). The magnitude of out-of-plane motion is approximately equal to the magnitude of extension predicted by the divergence in limb displacement directions. For the Bolivian orocline, recognized mechanisms for out-of-plane motion are regional paleomagnetic rotations and mapped strike-slip displacements along faults (Fig. 4). The objective is to restore both in-plane and out-of-plane deformation while minimizing the total kinematic incompatibility as quantified by the total overlap area at each time step. We test a range of permissible regional rotation magnitudes based on uncertainties in the available paleomagnetic and GPS data. Based on those results, additional strike-slip displacements are imposed along structures with geologic evidence for fault-parallel motion. The imposed displacements are kept at the minimum magnitude necessary to resolve large kinematic incompatibilities given the lack of direct constraints on strike-slip fault offsets. One limitation of the model is that polygon dimensions are fixed so the fault blocks are assumed to be rigid, ignoring internal strain. While this is a major simplifying assumption, recent work in the central Andes has shown that internal strain is negligible for the model region (Eichelberger and McQuarrie, 2015).

The fault-block boundaries are defined based on fault traces and fold axes from published geologic maps that range in scale from 1:100,000 to 1:500,000 (Fig. 4A) (e.g., Kley, 1996; Sheffels, 1988; McQuarrie, 2002; Müller et al., 2002; McQuarrie et al., 2008; Eichelberger et al., 2013). The blocks are regionally grouped into established physiographic zones present at the Bolivian orocline. From west to east, they are the Western Cordillera, Altiplano basin, Eastern Cordillera (EC), Interandean zone (IA), and Subandes (SA) (Fig. 1) (Isacks, 1988; Kley, 1996). The Western Cordillera is the modern volcanic arc with peak elevations of ∼6 km and forms the western boundary of the ∼3.8-km-high, internally drained Altiplano basin (Isacks, 1988). Due to extensive volcanic cover in the Western Cordillera, shortening is unconstrained and not included in the reconstructions presented here. The eastern limit of the Altiplano basin is defined by the EC, which has peak elevations of 5–6 km and is composed primarily of Paleozoic and Mesozoic rocks deformed by west-vergent folds and faults at the Altiplano basin–EC boundary. Structural vergence direction in the EC switches from west (back thrust) to east (fore thrust) toward the EC–IA border (Fig. 5). The IA is characterized by 2–3 km elevations and dominantly Ordovician–Devonian strata involved in thin-skinned deformation (Kley, 1996; McQuarrie, 2002; Müller et al., 2002) (Fig. 1B). The SA is marked by elevations of <2 km and an 8–14 km structural elevation decrease from the IA (McQuarrie, 2002; Eichelberger et al., 2013) (Fig. 5). Geophysical data indicate that the change in structural elevation from the SA to IA is accommodated by basement deformation below the IA (Kley, 1996; Schmitz and Kley, 1997; Wigger et al., 1994; Dunn et al., 1995; Zapata and Allmendinger, 2000). The SA is structurally characterized by broad, linear synclines that form valleys containing Tertiary synorogenic sediments (Figs. 1B and 5). The synclines are separated by ridge lines supported by faulted anticlines involving Devonian through Cretaceous strata (Baby et al., 1992; Dunn et al., 1995; McQuarrie, 2002). Elevations decrease through the SA into the Chaco foreland basin where topographic relief diminishes to a low foreland slope (Barnes and Heins, 2009). The boundaries of each zone (Fig. 1B) are based on changes in structural elevation, proposed involvement of basement thrusts at depth, and changes in fault vergence direction (see Kley, 1996; McQuarrie and DeCelles, 2002; McQuarrie, 2002). To allow for variable restoration vectors and deformation timing along the Cordillera, the model includes artificial boundaries perpendicular to structural trend that separate the model into five transects (Fig. 4A). The transect boundary positions are based on the map-view extent of regions with balanced cross-section shortening estimates, but have no structural significance.

The principal kinematic constraints for the Bolivian central Andes come from shortening estimates from balanced cross sections, shortening directions from mapped fault and fold orientations, deformation timing from thermochronology and basin studies, and vertical axis rotations documented in paleomagnetic and GPS data. In addition, geologic mapping and earthquake focal mechanisms indicate strike-slip displacements at the orocline axis along extensive fault zones (Dewey and Lamb, 1992; Funning et al., 2005), but direct constraints on the displacement magnitude are not available.

The magnitude of displacement for each block is assigned based on fault offsets as shown in the cross sections associated with that transect (Table 1; Fig. 4A). The southernmost transect, unlike those to the north, lacks a single continuous line of section across the fold-thrust belt (transect 5, Fig. 4), so the total displacement is based on a composite shortening estimate from structural studies in each of the physiographic units (SA—Dunn et al., 1995; IA—Kley, 1996; EC—Müller et al., 2002). For fault blocks in all six transects, the shortening is restored in sequence, from east (SA) to west (EC). During restoration, active physiographic units are internally restored, indicated by expanding space between fault blocks. Outboard units to the west that have earlier deformation periods (restored in subsequent time steps) are passively displaced with no intraunit separation between fault blocks. Once a unit has been fully restored, the constituent fault blocks remain fixed. The primary challenge in reconstructing the orocline is that the curvature of the fold-thrust belt implies divergent plane-strain displacement paths that intersect as shortening is restored. Shortening magnitude estimates used for the reconstruction are quite similar and range from 276 km (40%) in the distal northern limb (McQuarrie et al., 2008) to 327 km (37%) in the southern limb (Table 1) (McQuarrie, 2002). The northwest-southeast structural fabric of the northern orocline limb implies a northeast principal shortening direction while the southern limb is characterized by a north-south structural fabric and east-directed shortening (Table 1; Fig. 4A).

Timing of deformation in the model is based on a combination of thermochronologic cooling ages from each of the model transects and geologic constraints from mapping and basin sedimentation histories. The available thermochronology data varies by transect, it generally includes ⁴⁰Ar/³⁹Ar, zircon fission track, zircon (U-Th)/He, apatite fission track, and apatite (U-Th)/He ages from various stratigraphic depths and structural positions across all three major physiographic units (Fig. 5). These cooling ages have been used as a timing proxy for central Andes deformation under the assumption that cooling is driven by erosional exhumation in response to the relief generated by faulting (e.g., Barnes and Ehlers, 2009). Exhumation ages in each physiographic unit often overlap within uncertainty, preventing a fault-by-fault restoration sequence. The broad overlap in exhumation ages has been interpreted to reflect periods of distributed deformation both within each physiographic unit and across unit boundaries (Barnes et al., 2006, 2008; McQuarrie et al., 2008; Eichelberger et al., 2013). Rather than assuming a structure-by-structure propagation sequence, fault blocks within each physiographic unit are simultaneously restored over the same time period, beginning with the earliest exhumation age in that region (Fig. 5; Table 1). Overall, the oldest exhumation onset ages are found in the EC and generally become younger eastward into the SA (Barnes et al., 2006, 2008, 2012; McQuarrie et al., 2008; Eichelberger et al., 2013). While thermochronology can be used to interpret the onset of deformation within a physiographic region, determining the ages of cessation is more difficult. Geologic constraints show that deformation in the northern limb EC had ceased by 25 Ma (Gillis et al., 2006; Horton, 2005) while GPS data and geomorphology indicate that the SA is actively deforming today (Brooks et al., 2011). There are no data that directly constrain the termination of IA deformation. However, the kinematic systems proposed in balanced cross sections suggest that IA deformation would have ceased by the time SA deformation was underway, although IA exhumation would continue (McQuarrie et al., 2008) (Fig. 5). In general, deformation in the SA is restored from 15 to 0 Ma across the entire orocline. Timing of IA exhumation varies by location, spanning 30–5 Ma in the northern limb, 40–15 Ma at the orocline axis, and 25–10 Ma in the southern limb. As a result, IA deformation is restored in the model beginning with the earliest exhumation age and ending at 15 Ma, when SA deformation begins. In the EC, fault blocks are restored over 45–25 Ma in the northern limb, 50–20 Ma at the orocline axis, and 50–20 Ma for the southern limb.

At the regional scale, paleomagnetic and GPS data from Peru to Chile resolve the central Andean rotation pattern system of counterclockwise rotations in the northern limb and clockwise rotations for the southern limb (Figs. 3 and 4). At smaller scales, the rotations are subject to local structural influences, have large uncertainties, and were sampled from lithologies older than Late Cretaceous that may record rotations unrelated to Andean deformation. In addition, the data are concentrated in volcanic forearc rocks but sparse in the portion of the retroarc fold-thrust belt covered by the reconstruction. Collectively, these factors make it difficult to identify rotations related to the development of orogenic curvature. To determine representative magnitudes of rotation related to recent Andean deformation, we focus on paleomagnetic data from synorogenic sediments located in the SA (compiled by Roperch et al., 2006). The SA definitively accommodated the most recent phase of deformation in the central Andes from ca. 15 Ma to present. SA structures are dominantly fault-propagation folds with consistently linear local orientations, typical of foothills-style deformation (Dahlstrom, 1969). There is little evidence for structural complexities that could result in local block rotations, suggesting that rotations recorded by Tertiary synorogenic sediments document regional rotation. In addition, GPS data (Allmendinger et al., 2005) provide additional constraints on regional rotations. Extrapolating modern rotations over geologic time assumes constant rotation rates; however, at a regional scale, the convergence direction for the subducting Nazca plate has been constant since SA deformation began ca. 15 Ma (Pardo-Casas and Molnar, 1987; Somoza, 1998). Modeling the rotational component of oroclinal bending as contemporaneous with shortening in the SA differs from interpretations that attribute bending of the Andes to regional north-south compression unrelated to fold-thrust belt shortening (e.g., Johnston et al., 2013). Given that the SA in the model area has been actively exhuming since at least 5 Ma (Barnes et al., 2006, 2008, 2012; Eichelberger et al., 2013), continues to shorten today (Brooks et al., 2011), and modern backarc rotation rates are low (≤2° m.y.^–1; Allmendinger et al., 2005), it is unlikely that rotational bending postdates SA deformation.

Due to the large uncertainties in the paleomagnetic data, we perform two reconstructions to test minimum and maximum magnitudes of limb rotation. Paleomagnetic rotations in the southern limb SA between 21.5°S and 18°S range from –1.1° ± 14° (counterclockwise) to 11.1° ± 22° (clockwise) (Fig. 3). Two SA data points are available for the northern limb, recording 0.2° ± 21° at 15°S and –7.4° ± 7.2° at 12.7°S. These measurements are statistically equivalent, but within the context of the central Andean rotation pattern, we assume that –7.4° ± 7.2° is consistent with orocline scale deformation and complimentary to south limb rotation. The minimum rotation scenario is based on the mean SA rotation (∼7°; Table 1; Fig. 4B). The maximum rotation scenario is based on the maximum permissible rotation within uncertainty for the SA paleomagnetic data in the study region. In the reconstruction, the maximum rotation applied to the limbs is 13° (Table 1; Fig. 4B). The rotational range represented by the minimum and maximum rotation scenarios cover retroarc GPS rotations extrapolated over 15–0 Ma (south limb: 12°; north limb: –6°; Allmendinger et al., 2005) as well as the range of EC and Altiplano basin paleomagnetic rotations from the reconstruction region (Fig. 4C). Paleomagnetic rotations in the EC region of the orocline core are variable, likely due to the influence of possible strike-slip faulting in the area of Cochabamba where the available samples are located. Minimal rotations have been predicted for this region (Kley, 1999; Arriagada et al., 2008); however, the northwest-southeast–trending EC structures at the core are continuous with northern limb structures that undergo counterclockwise rotations. For the initial maximum and minimum limb rotation models, the core EC rotates with the northern limb while the IA and SA do not rotate (Fig. 4B).

The most prominent effect of restoring regional rotations is the progressive impact it has on changing the displacement directions defined in each region. The reconstruction assumes that rotation was a function of SA deformation: blocks that underwent less displacement also underwent a proportionally lower magnitude of rotation. Within the model framework, the magnitude of rotation restored at each time step during SA deformation is scaled by a given fault block’s displacement magnitude. Therefore, as with shortening, rotations are restored in sequence from the SA (recent deformation) back toward the EC (earlier deformation). Additional rotations can be added at earlier time steps to resolve local kinematic incompatibilities, but are not required because SA rotations account for measured EC rotations within uncertainty (Fig. 4C). By the time the SA is fully restored, blocks at the western edge of the SA and beyond undergo the full magnitude of assigned limb rotation. The incremental rotations at each time step then alter fault-block displacement paths in all subsequent time steps. The regional displacement direction in the minimum rotation model changes from 270° (west) to 262° (approximately west-southwest) in the southern limb and from 229° (approximately south-southwest) to 237° (approximately southwest) in the northern limb. The maximum rotation model displacement directions change to 256° (approximately west-southwest) in the southern limb and from 243° (approximately southwest) in the northern limb. Independently restoring the limb rotations and shortening translations in each transect imposes different poles of rotation for north and south limbs (similar to the external pivot model of Sussman et al., 2012). The poles are located away from the orocline axis and outside the reconstruction reference frame. Consequently, the only fixed points in the area of the reconstruction are along the modern foreland basin.

Of all the kinematic components documented at the Bolivian orocline, the magnitude of lateral displacements approximately parallel to fault trend is the least well constrained. Geologic and seismologic evidence shows that at least two faults have accommodated transpressional deformation in the region of the orocline axis and continue to do so today. The north-south–trending Rio Novillero fault (RNF, Fig. 6) at the orocline core was initially located based on the change from northeast-vergent thrust faults that characterize the north limb structural fabric to east-vergent structures characteristic of the south limb (Dewey and Lamb, 1992). This change occurs abruptly across several structures that have accommodated north-south–oriented, right-lateral motion indicated by locally offset lithologic contacts, but no direct constraints on the total north-south displacement are available (Eichelberger et al., 2013). Recent seismicity in the area also shows that the RNF accommodates right-lateral motion on nearly vertical faults planes today (Fig. 6) (Dewey and Lamb, 1992; Funning et al., 2005). The Cochabamba fault has been mapped as both a normal and left-lateral fault (Dewey and Lamb, 1992; Kennan et al., 1995; McQuarrie, 2002; Sheffels, 1995) due to the presence of two apparent pull-apart basins with bounding ranges that are 2–3 km higher in elevation (Cochabamba fault, CF, Fig. 6). Earthquake focal mechanisms are consistent with left-lateral fault displacements (Dewey and Lamb, 1992), but there are no geologic observations that constrain offset magnitude. Together the RNF and CF are proposed to form a conjugate fault system (Dewey and Lamb, 1992) that structurally separates the orocline core from the limbs, similar to the orocline model of Kley (1999). The extent of both faults beyond the orocline core region is ambiguous, but for different reasons. The regional change in structural trend associated with the RNF is traceable for ∼100 km until the structural trend becomes uniformly north-south in the southern limb (Fig. 6). The southern limb extension of RNF in the model follows east-vergent thrust faults where bedding orientations are slightly oblique to fault trend, similar to the structural relationship established to the north (Fig. 4A). The extent of the CF is best defined by the west-east–trending pull-apart basins, oblique to the northwest-southeast structural fabric. These features are only discernable over ∼50 km, suggesting that if the CF extends westward into the northern limb and Altiplano, fault-lateral offsets are probably distributed along the northwest-southeast–trending thrust faults that characterize the region west of the Cochabamba basins. In the model, the CF is extended westward following mapped thrust faults that are slightly oblique to the dominant trend but become parallel to the regional trend closer to the Altiplano (Fig. 4A). Where the RNF and CF parallel the local fold-thrust belt trend, any lateral displacement would be convolved with fault-normal shortening, making the lateral component of slip difficult to detect in the field. Lateral displacements are not mapped in the limbs, which suggests that the total offset on either fault system may be small and distributed over multiple faults.

In all of the reconstructions presented here, the kinematic incompatibility of a given time step is interpreted based on the total overlapping polygon area at that time step. The total overlap area at the final time step is the cumulative overlap from the entire model. However, for a physiographic region in a given transect, the maximum overlap magnitude may occur at an earlier time step. For example, the magnitudes of SA overlap in each model peaks when SA shortening has been fully restored then remains fixed in subsequent time steps. Since the IA and EC have concurrent restoration periods (ca. 45–20 Ma, depending on transect), peak overlap magnitudes may decrease as shortening is restored in adjacent zones and the area between polygons increases. In terms of improving the kinematic model by including fault-parallel displacements to minimize the overlap area, the location and time of overlap is as important as the magnitude. However, changes in a physiographic region’s overlap magnitude between time steps can reveal if imposing out-of-plane motions actually improves kinematic compatibility or if it simply transfers incompatibilities to other regions. To summarize the results of each reconstruction, we report the total overlap from the final time step as the total kinematic incompatibility and use it to compare reconstruction results. For regions with imposed fault-parallel displacements, local overlap magnitudes (overlap within a physiographic region of a transect) are compared at individual time steps to analyze whether out-of-plane movement resolves or redistributes kinematic incompatibility.

The magnitude of overlap area produced in a reconstruction using only plane-strain restoration constraints provides a benchmark for determining relative improvement in kinematic compatibility with increasing kinematic complexity (Fig. 7). Without rotation or fault-parallel displacements, the total kinematic incompatibility at the final stage of the restoration is 14,160 km². Of the total overlap, the majority (11,400 km², 80%) results from the intersection of the southwest-directed restoration path of the core with the west-directed south limb restoration path. No kinematic incompatibilities arose in the north limb. There are secondary regions of kinematic incompatibility in the orocline core that peak prior to the final time step due to the sequence of restoration. SA overlap (1000 km²) peaks at 15 Ma and the EC overlap (2500 km²) peaks at 25 Ma (Fig. 8A). These kinematically incompatible areas at the orocline core are due to internal changes in structural trend that are not found in either of the limbs (Fig. 7).

The minimum limb rotation (±7°) reconstruction (Fig. 9) results in substantially less kinematic incompatibility with the final magnitude of overlap dropping to 5300 km², >60% less overlap area than the plane-strain model (Fig. 8B). The regional reduction in overlap magnitude is largely due to the inclusion of limb rotation, which alters the limb displacement paths so that they are closer to parallel. Kinematic compatibility was further improved by imposing ∼90 km of south-directed, fault-parallel slip along the RNF over 50–0 Ma to eliminate overlap between the core and south limb EC in transect 4 (Fig. 9; Table 2). However at the local scale, propagating slip on the RNF into the southern limb increases EC overlap in transects 5 and 6 (Figs. 4 and 9). As a result the greatest total overlap (∼7000 km²) occurs at 25–20 Ma as the RNF propagates south, but this overlap subsequently decreases as the southern limb EC is restored (Fig. 8B). Overall, the reduction in overlap from the plane-strain model outweighs any local increases, but the introduction of localized kinematic incompatibilities in the south limb from RNF slip indicates that it may not be kinematically viable for the fault to accommodate ∼90 km of north-south slip. Instead, larger limb rotations and/or the inclusion of slip on the CF may be required to resolve the incompatibility of the core and southern limb displacement directions. The kinematic incompatibilities within the orocline core show no reduction in the IA and SA but are reduced in the EC (Fig. 8B). This is because the northwest-southeast–trending EC structures are allowed to rotate with the north limb from 15 to 0 Ma whereas the IA and SA do not rotate. IA and SA overlap at the core both peak at 15 Ma, then IA overlap decreases as deformation is restored from 20 Ma onward (Figs. 4B and 9).

Increasing the magnitude of limb rotation to the upper permissible limit of ∼13° further improves kinematic compatibility (Fig. 10) with a final overlap area of 3200 km², an ∼80% decrease over the final overlap in the plane-strain model (Fig. 8C). With the larger limb rotations, ∼45 km of north-south slip on the RNF was required to resolve kinematic incompatibilities between the core and southern limb, roughly half the magnitude required by the minimum limb rotation model (Table 2). Lower RNF displacement also introduced less kinematic incompatibility in the south limb, as shown by peak south limb EC overlap of 1300 km² at 20 Ma (Fig. 8C) compared to ∼3000 km² at 35 Ma in the minimum rotation model (Fig. 8B). Overall, the larger limb rotations greatly improved the kinematic compatibility of the core and southern limb while reducing the magnitude of slip required on the RNF, and by extension, the magnitude of overlap introduced in the southern limb. Within the orocline core, local kinematic incompatibilities related to changes in structural orientation remain unresolved in this model. The total overlap area for the maximum rotation model peaked at 15 Ma (∼4400 km²), at which time incompatibilities in the core IA and SA account for ∼60% (2560 km²) of the total overlap area (Fig. 8C). The next largest overlap component at 15 Ma is related to north-south slip on the RNF in the south limb EC (∼1200 km²) and composes 30% of the total 15 Ma overlap (Fig. 8C). Further improving the kinematic compatibility of the reconstruction will require a solution that can account for the changes in structural trend at the orocline core and further reduce the amount of RNF slip required to mitigate overlap between the southern limb and orocline core.

Overall, the maximum rotation model had the lowest magnitude of overlap at every time step (Fig. 8C) and imposed far lower strike-slip displacements on the RNF. Regional limb rotation resulted in displacement paths that were closer to parallel, resolving many of the regional incompatibilities related to intersecting plane-strain shortening directions (Fig. 11). In turn, this reduced the magnitude of fault-parallel displacement interpreted at the orocline core in order to resolve kinematic incompatibility between the core and southern limb. We believe that this indicates that at least ∼13° of rotation in both limbs is required to achieve greatest degree of kinematic compatibility in a reconstruction with the current constraints available. This is slightly larger than the 5°–10° of rotation estimated by Kley (1999), but agrees with GPS rotation magnitudes extrapolated over the 15–0 Ma SA deformation window (Allmendinger et al., 2005). For the region of the fold-thrust belt covered by this reconstruction (14.6°S–22.5°S), 13° limb rotations used are within the range of solutions proposed by Arriagada et al. (2008).

Even in the maximum rotation model, localized areas of overlap persist at the orocline axis despite the regional improvement in kinematic compatibility. In addition, introduced kinematic incompatibility related to slip on the RNF in the southern limb suggests that large right-lateral offsets may not be viable there. The CF was considered fixed for the minimum and maximum rotation models, but may have accommodated transpressional deformation based on its oblique orientation relative to the surrounding EC structural fabric (Fig. 6). Further refining the core kinematics to reduce the remaining overlaps beyond the available constraints can resolve local incompatibilities as well as elucidate the role of the CF in accommodating deformation at the core.

To investigate these possibilities, we present an additional, preferred model that builds on the results from the maximum rotation model (Fig. 12). For the orocline limbs, displacement magnitude, limb rotation, and timing are all kept the same. At the orocline core we prescribe SA rotation to be equivalent to limb rotation magnitudes and impose small rotations in the core IA and EC as well as fault-parallel offsets on the CF to resolve remaining internal core overlaps. As with the EC, northwest-southeast–trending SA fault blocks rotate with the northern limb while more north-south–trending SA blocks rotate with the southern limb (dashed rotation lines, Fig. 4A). The final overlap with these adjustments is ∼1600 km² with a peak overlap of 1900 km² at 15 Ma for a ∼50% reduction in peak and final overlap from the maximum rotation model (Fig. 13).

By allowing the core SA to rotate the full magnitude of limb rotation, the internal overlap decreases slightly from 900 km² to 700 km². Slight differences in the structural orientation of the SA in the core and southern limb result in a shift of some SA overlap to the southern limb, increasing overlap from <100 km² in the maximum rotation model to ∼600 km². The total magnitude of SA overlap remains similar (1400 km² compared to 1000 km² in the maximum rotation model), indicating that the internal SA overlap was redistributed southward (Fig. 13). The remaining SA kinematic incompatibility at the core is related to the gradual change in SA structural orientation in the area, implying local, nonparallel displacement directions (Fig. 4A). A possible solution would be small but distributed fault-parallel displacements along the north-south–trending SA structures, similar to the model of Kley (1999).

As with the first two models, northwest-southeast–trending EC blocks defined by structures west of the RNF rotated with northern limb but at a lower magnitude (∼5° instead of 13°; Fig. 4A). The EC and IA structures east of the RNF had no defined rotation in the first two models, but here are rotated ∼5° with the southern limb based on their close alignment with the north-south southern limb fabric. In effect, the RNF becomes a rotational boundary that accommodates transpression related to the counter rotation of the two structural domains. By allowing the north-south EC and IA structures at the core to counter rotate, we resolve all the internal core overlap related to the east to west change in structural orientation between the SA, IA, and EC (cf. IA overlap in Figs. 8C and 13). Slip along the CF was applied to accommodate both the counter rotation of the northwest-southeast and north-south structural domains and shortening within the IA at the orocline core. By allowing ∼50 km of left-lateral CF slip over 25–0 Ma (Table 2), kinematic incompatibility between the EC at the orocline core and the southern limb was largely eliminated. This is equivalent to the magnitude of right-lateral RNF slip from the maximum rotation model, but has the advantage of not producing additional overlap in the limbs. The majority of orogen-parallel strike-slip motion occurs on the CF in this model while the RNF accommodates <10 km of right-lateral motion.

In the preferred model the RNF and CF behave as a conjugate fault system (Dewey and Lamb, 1992). However, rather than accommodating purely strike-slip motion, the reconstruction suggests that both accommodated out-of-plane transpressional deformation that was not previously recognized (Fig. 12). For the RNF, the preferred model estimates ∼40–50 km of strike-normal shortening, which we have distributed across 4 fault blocks between the RNF and the IA-EC boundary over 25–15 Ma and concentrated along the RNF trace from 15 to 0 Ma (Table 2). The preferred model estimates ∼60–70 km of strike-normal shortening on the CF from 45 to 0 Ma (Table 2). This amount of shortening is a function of how the IA is restored at the orocline core from 45 to 15 Ma, but is currently not accounted for in shortening estimates in the region. Less than 20 km of CF shortening occurred from 15 to 0 Ma, because the IA and EC were inactive. As a result, motion on the CF was dominated by the strike-slip displacements imposed to resolve kinematic incompatibility. The geologic map patterns for both faults (McQuarrie, 2002; Eichelberger et al., 2013) do not rule out the possibility of significant shortening, but most thrust faults in the central Andes have offsets <15 km (e.g., McQuarrie, 2002). This suggests that transpression is most likely distributed over several structures rather than localized to the degree modeled here. The CF region has higher structural elevations and deeper exposed stratigraphy than the RNF region, consistent with localized shortening. The CF also features large transtensional basins with 2–3 km of relief, indicating nontrivial offsets along the fault. Comparatively, the RNF has only minor topographic expression and is largely defined by the change in structural orientation across the apparent fault trace.

The CF and RNF are both located within the EC, which is restored over 50–25 Ma in this reconstruction (Fig. 12). As such, the transpressional displacements applied to the faults over 20–0 Ma in the preferred model represent out-of sequence deformation relative to the period of deformation interpreted from thermochronology and geologic constraints. Apatite fission track cooling ages from the region indicate the time at which the rocks were exhumed to ∼110 ± 10 °C (Gallagher et al., 1998). As a result, later phases of deformation may not be resolved well, especially if the structures in question were not sampled in detail, or accompanied by significant exhumation. The geologic constraints that led us to terminate EC deformation (Gillis et al., 2006; Horton, 2005) do not strictly prohibit localized deformation in the vicinity of the RNF or CF. Modern seismicity indicates that both faults continue to accommodate strike-slip displacements on nearly vertical fault planes (Fig. 6), but the available focal mechanisms do not indicate active thrust faulting (Dewey and Lamb, 1992; Funning et al., 2005). From 5 to 0 Ma in the preferred model, the CF and RNF accommodate 5–10 km of fault-parallel slip, but negligible strike-normal shortening (<5 km; see 5 Ma time step in Figs. 10 and 12), consistent with modern seismicity.

Both the maximum rotation and preferred reconstructions demonstrate that out-of-plane displacements are critical to obtaining the most self-consistent reconstruction with the available kinematic constraints. Regional limb rotations and transpressional deformation focused at the orocline core have the combined effect of increasing shortening estimates (Pueyo et al., 2004) over previously published 2-D plane-strain estimates alone (e.g., Müller et al., 2002; McQuarrie, 2002). The preferred model estimates ∼340 km of shortening at the orocline axis, compared to ∼270 km estimated by cross section alone (Eichelberger et al., 2013). The ∼70 km of additional map-view shortening at the core is due to the combined effect of transpressional shortening (∼50 km) at the RNF and ∼20 km of additional SA shortening due to the imposed curved slip paths from limb rotation. Including Altiplano shortening (∼40 km; McQuarrie, 2002) and the shortening estimate error at the orocline axis (∼15%; Eichelberger et al., 2013), total map-view shortening at the orocline axis would be 380 ± 50 km. Total map-view shortening for both orocline limbs is also higher compared to published plane-strain estimates (Fig. 14). For the northern limb, the preferred model predicts ∼370 km compared to ∼280–300 km (McQuarrie, 2002; McQuarrie et al., 2008). There, the 70–90 km increase is due to ∼20 km of additional SA shortening imposed by rotation and ∼60–70 km of transpressional shortening along the CF (Table 2). In the southern limb the preferred model predicts ∼300–350 km of shortening compared to ∼280–330 km (Dunn et al., 1995; Kley, 1996; McQuarrie, 2002; Müller et al., 2002). The ∼20 km increase in southern limb shortening is related to SA rotation because RNF offset is negligible in the preferred model.

Shortening rate ranges predicted by the models are 5–10 mm yr^–1 in the EC, 2–4 mm yr^–1 in the northern and central IA, 6–10 mm yr^–1 in the southern IA, and 5–9 mm yr^–1 in the SA (Table 3). The predicted shortening rate during EC deformation (ca. 50–20 Ma) are consistent with rates calculated from balanced cross sections (Barnes et al., 2008; McQuarrie et al., 2008) but are at the upper end the 0–8 mm yr^–1 rate proposed for early ca. 45–30 Ma EC deformation in southernmost Bolivia (Oncken et al., 2006) (transect 5, Fig. 4A). SA shortening rates from 15 to 0 Ma in the models are also similar to cross-section shortening rates, but lower than the 9–16 mm yr^–1 SA rates from 10 to 0 Ma proposed by Oncken et al. (2006). Although the modeled shortening rates do not consider short-term (<5 m.y.) fluctuations, the modeled 5–9 mm yr^–1 SA rate is more consistent with Quaternary shortening rates of 7–10 mm yr^–1 (Echavarria et al., 2003; Uba et al., 2009). The maximum rotation and preferred models impose limb rotations similar to GPS data (Allmendinger et al., 2005) and predict average SA shortening rates of 6–9 mm yr^–1, at the low end of 9–13 mm yr^–1 GPS rates from the southern SA (Brooks et al., 2011).

In the minimum rotation model SA shortening is 5–10 km greater than cross-section estimates and 20–25 km greater in the maximum rotation and preferred models (Table 2). At face value, the preferred model estimates increase the SA shortening differential between southern Peru to the orocline axis to 88 km (based on 17 km of SA shortening; Gotberg et al., 2010), corresponding to ∼6° rotation for the northern limb. Between the axis and northern Argentina, the predicted differential is 45 km (60 km of SA shortening; Echavarria et al., 2003), corresponding to ∼5° of rotation. In both cases, this represents slightly less than half the total rotation imposed in the model (∼13°). If curved SA slip paths accommodated limb rotation as modeled, regional bending unassociated with shortening may be partially responsible for modern orogenic curvature. Alternatively, if SA cross-section shortening estimates represent the total SA displacement field, only ∼3°–5° of rotation can be attributed to an SA shortening gradient (Eichelberger et al., 2013), while the remaining 8°–10° would be related to regional bending. Ultimately, the models indicate that ∼13° of limb rotation produces the most kinematically viable reconstructions, but the relative rotational contributions of differential shortening, curved slip paths, and regional bending are uncertain.

At the orocline axis, predicted map-view shortening from the preferred model is equivalent, within error, to the 400 km shortening estimate modeled by restoring the central Andean rotation pattern paleomagnetic data (Arriagada et al., 2008) (Fig. 14). Both rotationally modeled and plane-strain shortening estimates have been used to estimate crustal thickness based on cross-section area (e.g., Gotberg et al., 2010; Kley and Monaldi, 1998; McQuarrie, 2002). The reconstructions presented here suggest that the out-of-plane contribution to shortening may be as high as ∼90 km. The map-view shortening estimates for both limbs and the orocline axis exceed the magnitude required to account for the modern cross-sectional area (Fig. 14) without requiring significant additional shortening in the forearc. As a result, shortening at the orocline may be sufficient to account for, and even exceed, modern crustal thickness. From a crustal budget perspective, this opens the possibility that crustal material may have been transferred by lower crustal flow or removed entirely via delamination events. Enhanced lithospheric thickening related to orocline formation has been linked to delamination in the Cantabrian orocline of northern Iberia (Gutiérrez-Alonso et al., 2004).

The map-view reconstructions presented here integrate plane-strain shortening estimates, structural orientations, paleomagnetic rotations, and deformation timing constraints from thermochronology and geology into 3-D kinematic models for the Bolivian orocline. These models indicate that material displacements parallel to orogenic trend are critical to produce kinematically viable orocline reconstructions. In particular, out-of-plane displacements are largely the result of convergent orocline limb rotations that are accommodated at the orocline core by strike-slip and transpressional faulting. Overall, the reconstructed displacement field suggests that map-view shortening estimates may account for, or even exceed, modern crustal thickness at the central Andes. In general, plane-strain studies of deformation (such as balanced cross sections) in areas of orogenic curvature only record a portion of the total orogenic displacement field. Here, map-view shortening estimates from the reconstruction may exceed those required to account for modern crustal thicknesses. If correct, this suggests that formation of the Bolivian orocline contributed to localized crustal thickening and lower crustal loss proposed to explain rapid surface uplift in the Altiplano.

Our main conclusions are as follows.

1. By accounting for limb rotation, kinematic compatibility as measured by map-view overlap was drastically improved over a reconstruction with no rotation (plane-strain model). North-south displacements on the strike-slip RNF further reduced kinematic incompatibilities at the orocline core. Final overlap in the plane-strain reconstruction was ∼14,000 km² compared to ∼5000 km² in the minimum rotation model (6°), and ∼3000 km² in the maximum rotation model (13°). The minimum rotation model predicted ∼90 km of strike-slip displacement on the RNF fault while the maximum rotation model predicted ∼45 km of slip.

2. We developed a third, preferred model based on the maximum rotation model. The preferred model includes minor rotations at the orocline core and slip on both the RNF and CF, resulting in ∼1600 km² of final overlap. This model applied ∼50 km of left-lateral slip on the CF and <10 km of slip on the RNF. Imposed transpressional deformation predicted ∼60–70 km of additional strike-normal shortening in the vicinity of the CF and ∼50 km of shortening near the RNF. The model predicts that some of this shortening occurs after the main phase of deformation in the EC (before 25 Ma) where the faults are located.

3. In the preferred model, the combined effect of limb rotation and transpressional deformation produces map-view shortening estimates that are greater than in cross-section shortening estimates. Map-view shortening predicted by the preferred model is ∼370 km in the northern limb, ∼380 km at the orocline core, and 320–340 km in the southern limb. Lower shortening estimates are predicted for the southern limb because transpressional shortening is limited to the orocline core and northern limb. The 70–90 km increase in shortening in the northern limb and orocline core is the result of 60–70 km of transpressional shortening and ∼20 km of additional SA shortening due to curved slip paths.

4. Increases in SA shortening due to the curved slip paths imposed in the model would result in 5°–6° of limb rotation due to differential shortening between southern Peru, the orocline core, and northern Argentina. This is slightly less than half the total rotation applied in the preferred model, suggesting that regional bending may be a factor. If limb rotation is not the result of curved slip as modeled here, 8°–10° of the total 13° limb rotations required for kinematic compatibility would have to be the result of regional bending.

5. Map-view shortening at the orocline axis is predicted to be 380 ± 50 km. This value exceeds the magnitude required to account for modern cross-section area of the crust. This suggests that deformation at the orocline may be sufficient to account for the modern crustal thickness without significant additional shortening in the forearc. This opens the possibility that any excess may have been transferred along strike by lower crustal flow or removed from the system entirely by processes such as delamination.

We thank Stephen Johnston and an anonymous reviewer for their feedback and insightful comments on the initial version of this paper. Early versions of the reconstructions benefitted from discussions with Carmala Garzione, Laura Wagner, Todd Ehlers, Susan Beck, George Zandt, Brian Horton, Chris Poulsen, and their respective students. These interactions took place at workshops organized as part of the Central Andean Geodynamics and High Topography (CAUGHT) project funded by the National Science Foundation (NSF grant EAR-0908972 to McQuarrie). Blair Schoene provided comments on an early version of the manuscript. Melissa Brenneman created the original ArcGIS (geographic information system) script that produced the reconstructions shown here.