Reconciling Spatial and Temporal Patterns of Cenozoic Shortening, Exhumation, and Subsidence in the Southern Bolivian Andes

The Bolivian Andes are an archetypal convergent margin orogen with a paired fold-thrust belt and foreland basin. Existing chronostratigraphic constraints highlight a discrepancy between unroofing of the Eastern Cordillera and Interandean Zone fold-thrust systems since 40 Ma and the onset of rapid sediment accumulation in the Subandean Chaco foreland after 11 Ma, previously attributed to Miocene climate shifts. New results from magnetostratigraphic, backstripping, erosional volumetric calculations, and flexural modeling efforts are integrated with existing structural and thermochronologic datasets to investigate the linkages between shortening, exhumation, and subsidence. Magnetostratigraphic and backstripping results determine tectonic subsidence in the Chaco foreland basin, which informs flexural models that evaluate topographic load and lithospheric parameters. These models show that Chaco foreland subsidence is consistent with a range of loading scenarios. Eroded volumes from the fold-thrust belt were sufficient to fill the Chaco foreland basin, further supporting the linkage between sediment source and sink. Erosional beveling of the Eastern Cordillera, local intermontane sediment accumulation after 30–25 Ma, and regional development of the high-elevation San Juan del Oro geomorphic surface from 25 to 10 Ma suggest that the western Eastern Cordillera did not store the large sediment volume expected from erosion of the fold-thrust belt, which arrived in the Subandean Zone after 11 Ma. Eocene to middle Miocene foreland basin accumulation was likely focused between the Eastern Cordillera and Interandean Zone, and has been almost completely recycled into the modern Subandean foreland basin. The delay between initial fold-thrust belt exhumation (early Cenozoic) and rapid Subandean subsidence (late Cenozoic) highlights the interplay between protracted shortening, underthrusting, and foreland basin recycling. Only with sufficient crustal shortening, accommodated by eastward advance of the fold-thrust belt and attendant underthrusting of Brazilian Shield lithosphere beneath the Subandes, did the Subandean zone enter proximal foreland basin deposystems after ca. 11 Ma. Prior to the late Miocene, the precursor flexural basin was situated westward and not wide enough to incorporate the distal Subandean Zone. These results highlight the interplay between a range of crustal and surface processes linked to tectonics and Miocene climate shifts on the evolution of the southern Bolivian Andes.

Despite the wealth of information on the tectonic development of the central Andes, a full integration of the varied structural, thermochronological, stratigraphic, and geomorphologic datasets remains challenging. Specifically, the lack of Cenozoic mass balance constraints highlights challenges in understanding the linkages among shortening, erosion, and deposition within this paired fold-thrust belt and foreland basin system (Barnes and Heins, 2009). Key discrepancies underscore this persistent problem. For example, although the Eastern Cordillera recorded protracted shortening spanning the Eocene to Miocene (McQuarrie, 2002;Barnes et al., 2006Barnes et al., , 2008Gillis et al., 2006;Ege et al., 2007;Anderson et al., 2018), Subandean basin fill recorded punctuated sediment accumulation since the middle to late Miocene, possibly attributable to climate shifts . This contrast highlights a mismatch between active deformational phases in the Eastern Cordillera fold-thrust belt and the main phase of subsidence in the accompanying foreland basin (Echavarria et al., 2003;Uba et al., 2006Uba et al., , 2009Calle et al., 2018).
The present study is motivated by key datasets and geologic relationships that facilitate investigation of shortening, erosional exhumation, and sediment accumulation. A series of accurate, incrementally balanced cross sections (McQuarrie, 2002;Anderson et al., 2017Anderson et al., , 2018 are complemented by rich thermochronologic records (Barnes et al., 2006Gillis et al., 2006;Ege et al., 2007;McQuarrie et al., 2008;Rak et al., 2017;Anderson et al., 2018) (Figure 1) that document the timing and location of contraction, and provide constraints on the magnitude of erosional unroofing. Relict foreland basin fill preserved within the orogenic interior (Kennan et al., 1995;Horton, 1998Horton, , 2005Mosolf et al., 2011), hinterland plateau (Sempere et al., 1990;Baby et al., 1997;Horton et al., 2001Horton et al., , 2002Murray et al., 2010), and Subandean Zone (Uba et al., 2006(Uba et al., , 2009Calle et al., 2018) determine the timing and magnitude of sediment accumulation. The San Juan del Oro surface (Figure 1) is an ancient geomorphic surface that remains relatively undeformed and is currently situated at ca. 3 km elevation (Gubbels et al., 1993;Hérail et al., 1996;Kennan et al., 1997;Müller et al., 2002). The surface marks the cessation of shortening across the Eastern Cordillera, and presents a datum against which to gauge erosion and incision since final surface construction at roughly 10 Ma. The colocation of these features provide a unique opportunity to investigate mass balance in the context of regional shortening and basin subsidence.
Here, multiple new datasets are combined: magnetostratigraphic results from a Miocene stratigraphic section in the Subandean Zone to facilitate a new backstripping analysis; estimates of the volume of material eroded from the Eastern Cordillera and Interandean Zone fold-thrust belt; and estimates of the eroded volume beneath the San Juan del Oro surface and the volume of sediment preserved in the Chaco foreland basin (as restricted to the southern Bolivian transect considered here). Coupled with existing results on the timing, magnitude, and kinematics of shortening in the Eastern Cordillera, Interandean Zone, and Subandean Zone, the following questions are addressed for southern Bolivia. (1) What are the timing and magnitude of tectonic subsidence in the Subandean foreland basin? (2) Is the magnitude of tectonic subsidence comparable to flexural modeling predictions?
(3) Are eroded rock volumes from potential source regions consistent with the volume of preserved foreland basin fill? (4) Why did the Subandean foreland basin experience rapid subsidence 10s of Myrs after the onset of fold-thrust belt contraction?

GEOLOGIC BACKGROUND
The Central Andean fold-thrust belt can be divided into distinct structural zones defined by major changes in topographic elevation, stratigraphic exposure level, and structural style (Kley, 1996;McQuarrie, 2002;Anderson et al., 2017). From west to east, these include the Eastern Cordillera, the Interandean Zone, the Subandean Zone, and the Chaco Plain (Figure 1). The Eastern Cordillera is divided into a backthrust belt in the west and a forethrust belt in the east, centered about a north-south axis through the Tupiza region (Kley et al., 1997;McQuarrie, 2002;Müller et al., 2002;Sempere et al., 2002;Anderson et al., 2017). The western backthrust belt is characterized by closely spaced faults and short wavelength folds (<2 km), along with fault-bounded late Oligocene to Miocene intermontane basins (Horton, 1998;Müller et al., 2002;Anderson et al., 2017). The eastern forethrust belt is characterized by more widely spaced structures, and the presence of two significant regional features: the Camargo syncline, and the ca. 55 km wide Cuesta de Sama anticlinorium. The latter persists for hundreds of kilometers along strike from northern Argentina to central Bolivia (Kley, 1996;McQuarrie, 2002;Anderson et al., 2017). East of the Cuesta de Sama anticlinorium, shallower structural levels are exposed in the Interandean Zone as the basal decollement cuts up from Cambrian-Ordivician rocks to a regional decollement within Silurian shales (Baby et al., 1992;Dunn et al., 1995). Similar to the Eastern Cordillera, the Interandean Zone can be divided into a western backthrust zone and eastern forethrust zone (Anderson et al., 2017). However, the structural style in the Interandean Zone is characterized by more closely spaced thrust faults accompanied by tight fault bend folds (1-3 km), stratigraphic exposures of predominantly Silurian and Devonian rocks, and a distinct lack of any Cenozoic synorogenic sediments (Kley, 1996;Anderson et al., 2017). The Subandean Zone is defined by a valley-and-ridge topography, with broad wavelength folds (5-10 km) bound on their eastern limbs by east-vergent thrust faults that are continuous along strike for hundreds of kilometers (Baby et al., 1992;Dunn et al., 1995;Anderson et al., 2017). Stratigraphic exposure levels transition from Silurian-Devonian rocks in the Interandean Zone to mainly Carboniferous-Neogene rocks in the Subandean Zone, with Carboniferous strata exposed in cores of anticlinal ranges and Neogene synorogenic rocks in the intervening synclinal valleys (Anderson et al., 2017). In the western Chaco plain, the active frontal thrust (Mandeyapecua thrust) lies ca. 60 km east of the Subandean topographic front (Figure 1) (Costa et al., 2006;Brooks et al., 2011). The structure of the intervening region is characterized by gentle folds and low-offset thrust faults concealed beneath low-relief Quaternary sedimentary cover (Baby et al., 1992;Dunn et al., 1995;Horton and DeCelles, 1997;Anderson et al., 2017).
Shortening in the Eastern Cordillera began in the Eocene. The oldest zircon (U-Th)/He cooling ages suggest exhumation and cooling of the Cuesta de Sama anticlinorium potentially as early as ca. 40-39 Ma. However, these samples have relatively large errors . Alternatively, thermal modeling of paired zircon (U-Th)/He, apatite fission track, and apatite (U-Th)/He samples indicate rapid exhumation of the Cuesta de Sama anticlinorium may have occurred as late as ca. 32 Ma (Ege et al., 2007;Anderson et al., 2018). Despite uncertainty regarding the precise timing and position of initial contraction, upper crustal shortening advanced westward across the Eastern Cordillera backthrust belt primarily in late Eocene and Oligocene time (McQuarrie and DeCelles, 2001;DeCelles and Horton, 2003;Ege et al., 2007;Anderson et al., 2018), and reached the San Vicente thrust along the Eastern Cordillera-Altiplano boundary by ca. 27 Ma (Müller et al., 2002;Ege et al., 2007;Anderson et al., 2018). The end of upper crustal shortening in the Eastern Cordillera is marked by diminished sedimentation rates and deposition of intermontane basin fill unconformably on deformed Paleozoic strata of the Eastern Cordillera (Horton, 2005). At this time, rapid cooling propagated eastward into the Interandean Zone from 25 to 18.5 Ma (Ege et al., 2007;Anderson et al., 2018). Depositional ages of Neogene strata determined from interbedded tuffs and apatite (U-Th)/He cooling ages show the thrust front migrated into the western Subandean Zone by ca. 8 Ma, propagated eastward to the Villamontes section by 1.5-2 Ma, and reached the active thrust front (Mandeyapecua thrust) after 2.1 Ma (Costa et al., 2006;Uba et al., 2009;Hulka and Heubeck, 2010;Lease et al., 2016;Anderson et al., 2018;Calle et al., 2018).
This study reports new magnetostratigraphic and backstripping results for the Angosto del Pilcomayo section (Uba et al., 2005) ca. 4 km west of Villamontes, Bolivia in the Chaco foreland basin (Figure 1). Along the banks of the Río Pilcomayo, ca. 1,370 m of exposed stratigraphic section spans the Petaca, Yecua, Tariquia, and Guandacay Formations. These Neogene sedimentary rocks represent retroarc foreland basin fill of the Chaco basin. The western third of the basin has been incorporated into fold-thrust structures of the Subandean Zone. The basin continues to the east where it onlaps distally onto the Brazilian shield. Sediment accumulation began in the late Oligocene synchronous with shortening in the Eastern Cordillera (Uba et al., 2005).
The synorogenic sedimentary wedge disconformably overlies Cretaceous eolian sandstones containing large-scale trough-cross bedding. The Angosto del Pilcomayo section is composed of the non-marine Petaca, Yecua, Tariquia, Guandacay, and Emborozú Formations (Figure 2) (Uba et al., 2005). A brief synthesis of these formations from Uba et al. (2006) follows. The Oligocene to middle Miocene Petaca Formation comprises fluvial sandstones, abundant pedogenic horizons, and reworked pedogenic conglomerates. The Miocene Yecua Formation has been interpreted as spanning marginal marine, lacustrine to fluvial conditions across the Subandean Zone. In the study area, the middle Miocene Yecua Formation overlies the Petaca Formation, possibly with a subtle low-angle unconformity, and hosts fluvial overbank and channel facies. The upper Miocene Tariquia Formation contains interbedded thin-and thickbedded sandstone, mudstone, and fine-grained sandstone, as well as poorly developed pedogenic horizons and abundant mudcracks. This has been interpreted as a product of a distal fluvial megafan environment. The upper Miocene to Pliocene Guandacay Formation contains conglomerate, sandstone, and mudstone consistent with braided streams in a medial fluvial megafan environment. The upper Pliocene to Pleistocene Emborozú Formation contains upward-coarsening thick-bedded conglomeratic strata displaying laterally continuous sheetlike bedding geometries attributed to deposition in a proximal fluvial environment.

Magnetostratigraphy
Presented here is an overview of the magnetostratigraphic sampling and analytical methods. A detailed summary is provided in the Supplementary Materials. In the field, individual beds were selected as potential sample sites on the basis of lithology and stratigraphic level. The stratigraphic position of sampled beds are anchored to dated volcanic tuff horizons (Uba et al., 2005(Uba et al., , 2006(Uba et al., , 2007(Uba et al., , 2009, providing foundational geochronologic constraints. Sampling focused on the Yecua, Tariquia, and lower ca. 100 m of the Guandacay Formation, which contain continuous exposures of siltstone to fine-grained sandstone lithologies best suited for magnetostratigraphic analysis. Samples were not collected from coarse-grained, discontinuous exposures intervals, including the underlying Petaca Formation and the overlying main body of the Guandacay or Emborozú Formations. Sample cores were obtained using a modified chainsaw with a circular 1-inch diameter, diamond tipped drill bit. At each sample site, a minimum of three separate cores (A, B, and C) were collected over a restricted stratigraphic range (generally less than 10-20 cm) and lateral distance (less than 1 m). The orientation of each core was measured before extraction using an orientation stage and compass. A total of 207 sites were collected from different stratigraphic horizons within the 1,370 m section, yielding an average stratigraphic spacing of ca. 6.5 m between successive sites. Using a customized trim saw with three circular, parallel, diamond tipped blades, each cylindrical core was cut into one or more half-inch thick disk-shaped specimens. To measure the remanent magnetization of the sample cores, each specimen was placed on the automated sample handling tray linked to a 2-G rock magnetometer in the Paleomagnetics Laboratory at The University of Texas at Austin. After measurement of the natural remanent magnetization (NRM), each specimen underwent a series of systematic thermal demagnetization steps using a shielded oven, measuring the thermal remanent magnetization (TRM) after each heating experiment. Due to instrument issues, only the A cores from all sites, and sites 1-136 of B cores could be analyzed. After progressive thermal demagnetization, the primary magnetization orientation for each specimen was analyzed using Paleomag X software. The resulting sequence of FIGURE 2 | Synthesis of magnetostratigraphic results. On the left, stratigraphic section geochronologic results from U-Pb geochronology at four positions (Uba et al., 2007) in the Angosto del Pilcomayo section near Villamontes, Bolivia. Site results from magnetostratigraphic investigation, showing the tilt-corrected declination of the primary magnetization vector determined for specimens at each site. N denotes normal polarity with a declination of 0 • , and R denotes reversed polarity with a declination of 180 • . Two possible interpreted polarity results are correlated with the Global Polarity Time Scale, pinned by the geochronologic ages shown. The ages and stratigraphic levels of the chrons in each correlation are plotted to define sediment accumulation curves. Correlation 1 is shown with white symbols, and average sediment accumulation rates are shown as linear interpolations (undecompacted). Correlation 2 is shown with light blue symbols. Two slightly different average sediment accumulation rates are shown, resulting in correlations 2.1 and 2.2. normal and reversed polarity intervals was then correlated to the Geomagnetic Polarity Time Scale (GPTS) (Gradstein et al., 2004), tied to stratigraphic horizons with known chronostratigraphic ages (Figure 2) (Uba et al., 2007). Two slightly different correlations (1 and 2) to the GPTS are presented. For correlation 2, average rates of deposition may be interpreted in two ways, leading to correlations 2.1 and 2.2.

Backstripping
The new magnetostratigraphic results were used to conduct 1D Airy backstripping with exponential reduction of porosity to determine the magnitude of tectonic subsidence at the Villamontes section. OSXBackstrip software by N. Cardozo was used to implement methods described in Watts (2001) and Allen and Allen (2013). The sections were modeled as nonmarine basins. Two scenarios were conducted, one with pure sandstone lithology, and one with pure shale lithology. This was done to investigate the relative importance of lithology on final magnitudes of tectonic subsidence, which are minimal (see section "Results" below). For sandstone, a density of 2,600 kg/m 3 , porosity coefficient c of 0.270 km −1 , and original porosity of 49% were used. For shale, density of 2,500 kg/m 3 , porosity coefficient c of 0.510 km −1 , and original porosity of 63% were used. Published ages and maximum stratigraphic thicknesses were used to model the Emborozú section . Tectonic subsidence was also estimated at two other positions east of the Villamontes succession along the line of section. These two additional locations lack chronostratigraphic control, and were used strictly to assess lateral variations in the total magnitude of tectonic subsidence. At these locations, we applied the ratio of tectonic subsidence magnitude to total sediment accumulation determined from the Villamontes section to determine equivalent tectonic subsidence at points B and C (Figure 3). The magnitude of tectonic subsidence is used for all three localities (Villamontes and point B and C) to constrain flexural modeling results.

Eroded and Deposited Volume
In Eastern Cordillera, Interandean Zone, and Chaco Foreland AreaErrorProp software (Judge and Allmendinger, 2011;Allmendinger and Judge, 2013) was used to calculate the crosssectional area of material eroded from the Eastern Cordillera and Interandean Zone, and the cross-sectional area of basin fill in the Chaco foreland. Eroded area was determined from the portions of the balanced cross section of Anderson et al. (2017) projected above the present surface. Eroded areas that belong to the Eastern Cordillera or the Interandean Zone are distinguished from each other and calculated separately based on the kinematic designation of each morphotectonic zone, as defined by Anderson et al. (2017Anderson et al. ( , 2018. In the Chaco foreland, the pre-Cenozoic basement slope was projected eastward to the approximate pinchout of Cenozoic basin fill in the distal Chaco foreland, constrained from the regional dip observed in seismic data and depicted in cross sections (e.g., Dunn et al., 1995;Tankard et al., 1995;Uba et al., 2009;Anderson et al., 2017). This projection agrees with maps of the foreland basin geometry (Horton and DeCelles, 1997). Once determined, cross-sectional areas were multiplied by 146 km, the strikeparallel (north-south) extent of the geologic map from Anderson et al. (2018), to determine the volume of material eroded from the Eastern Cordillera and Interandean Zone, as well as the volume of material preserved in the Chaco foreland basin. These calculations are restricted to the study area and do not apply to other latitudes of the orogenic belt. The 3D rugose nature of the basement-cover interface was not considered in the volumetric calculations of foreland basin fill. A north-south transect through the Subandean Zone based on subsurface datasets shows 1-2 km amplitude irregularities of this contact and an average depth of ca. 2.5 km depth below sea level (McGroder et al., 2015), although there is no systematic dip of the basement contact within the study area.

From Beneath the San Juan del Oro Surface
The geologic map of Anderson et al. (2018) was georeferenced in ArcMap to guide the interpolation of the high-elevation San Juan del Oro surface. The mapped distribution of the San Juan del Oro surface was extrapolated onto the MERIT DEM (multi-errorremoved improved terrain digital elevation model; Yamazaki et al., 2017). A slope map calculated from the MERIT DEM was used to aid mapping, which shows low-relief portions of the surface. By defining multiple points delimiting the preserved boundaries of the San Juan del Oro surface, the elevation was extracted from the MERIT DEM and added as an attribute field to the point feature class that defined the surface. Using these points, kriging interpolation was used to define the original low-relief surface geometry. The parameters for the kriging tool include the point feature class that indicates the boundary of the San Juan del Oro surface, and the z-value field, which indicates the attribute field in the point feature class that contains the elevation data. Additionally, the search radius, which defines the points used to interpolate the value for each cell of the output raster, was set at 100 points. All other kriging parameters were left as default values. Finally, to calculate the volume of material eroded from beneath the San Juan del Oro surface since its formation, the Raster Calculator tool was used to determine the difference between the interpolated San Juan del Oro surface raster and the MERIT DEM. All locations where the MERIT DEM was at a higher elevation than the kriged surface were excluded from these calculations.

Flexural Modeling
Matlab-based software FlexMC (Saylor et al., 2018) was used to perform Monte Carlo simulations that investigate a range of possible flexural loading configurations (variable dimensions and densities) and lithospheric conditions that could explain the magnitude of tectonic subsidence at three positions in the Bolivian foreland basin: Villamontes succession, and points B and C. The selected model space was defined as 800 km wide by 20 km high (10 km above and below a horizontal reference surface). The loads were allowed to vary between 50 and 300 km in width, and between 1 and 8 km in height. Elastic thickness could vary between 10 and 100 km. Three possible lateral patterns in effective elastic thickness (EET) were explored: laterally constant, linearly laterally variable, and exponentially laterally variable. Two possible geometries of the topographic load were explored: a rectangular load and a foreland-sloping load. The front (right) edge of the load was pinned at the 300 km position, and the load was situated to the left of the pin. The horizontal position of the three basin localities (Villamontes, B, and C) were measured from the Andean deformation front in the sequentially restored cross section of Anderson et al. (2018) at 8 Ma (Figure 3), and are depicted in the model (Figure 4; yellow points). Subsidence uncertainties of 500 m were implemented for the three basin localities. For each of the 6 model runs, 1,000 trials were conducted.

Magnetostratigraphy
Measured tilt-corrected declinations were used to determine magnetic polarity for 207 sites for all A specimens and sites 1-136 for all B specimens within the 1,370 m stratigraphic section (Figure 2). Using the magnetic polarity determined here and available geochronologic results from dated tuffs, two slightly different correlations 1 and 2 are presented. We define 30 polarity zones (15 normal, 15 reverse) in correlation 1, and 22 polarity zones (11 normal, 11 reverse) in correlation 2. These two correlations are consistent with independent U-Pb geochronologic results. Differences between correlations 1 and 2 come from the decision to include or exclude those polarity zones that are defined by a single site. In instances where a polarity zone is defined by a single site, this could be considered provisional until further tested. The two correlations are plotted to determine the sediment accumulation curve. Average sediment accumulation rates for intervals between inflection points in the curve are shown (Figure 2). For correlation 2, we explored the impact of selecting different inflection points (for example, between 800-1000 m) on average sediment accumulation rates. Despite this range of possible correlations and rates, there is general consistency among the different interpretations. The lowest interval (0-77 m) from both correlations shows an average accumulation rate of ca. 49 to ca. 77 m/Myr. The following interval, up to ca. 330 m, has an average accumulation rate of ca. 114 to ca. 140 m/Myr. The covered interval occupies the ca. 350-560 m level of stratigraphic section and prevents detailed magnetostratigraphic analysis of this part of the succession. Above the covered interval, the average accumulation rate ca. 314, ca. 677, and ca. 459 m/Myr for correlations 1, 2.1, and 2.2, respectively. This represents a 2 × (correlation 1) to 4.8 × (correlation 2.1) increase in sedimentation rate. Correlation 2.1 shows an average accumulation rate of ca. 677 m/Myr at ca. 8-7 Ma, consistent with the 628 m/Myr rate reported by Uba et al. (2007) for strata of the same age. The uppermost section, encompassing the 1,250-1,370 m interval, has an average accumulation rate of ca. 147 to ca. 233 m/Myr, which represents a slight reduction but is still higher than the rates for the lower interval. Overall, the sediment accumulation curve defined by new magnetostratigraphic data agrees with existing chronostratigraphic results from Uba et al. (2007), and provides higher resolution control on intervals between horizons with dated volcanic tuffs.

Backstripping
The magnetostratigraphic results define a series of magnetic chrons for the Villamontes section. The ages for the base and top of each chron, and associated stratigraphic level, were used for backstripping. The backstripping results show that the lithology chosen (sandstone or shale) predicts subtle differences in the decompacted sediment accumulation curves ( Figure 5). However, the difference in total tectonic subsidence between the two is minimal: 1,741 m (sandstone), and 1,759 m (shale). For the purposes of flexural modeling, a value of 1,750 m was used for the tectonic subsidence at the Villamontes succession. At the Emborozú section, results show sediment accumulation between 10.5 and 8.3 Ma only in the Tariquia and Emborozú Formations. The accumulation rate is similar to that observed in Pliocene-Pleistocene strata at the Villamontes section.
These results confirm the overall convex-upward sediment accumulation curve at the Villamontes section (Uba et al., 2007), consistent with flexural foreland basin deposition. At the incomplete Emborzú section, only a record of rapid sediment accumulation is preserved. At the Villamontes section, there is an inflection in the sediment accumulation curve at ca. 11-8 Ma, marking the transition to more rapid sediment accumulation in the foredeep depozone.

Eroded Volume
Eroded volumes (Figure 6) were calculated for the study area using the geologic map of Anderson et al. (2018). No attempt was made to determine whether the calculated volumes are appropriate for other segments of the Andes along strike to the north or south. No attempt was made to account for changes in volume and density associated with breakdown of solid rock into detritus, or the impact of compaction on volumetric changes after deposition. A comparison of density differences between sedimentary basin fill (2,500-2,600 kg/m 3 ) and potential Eastern Cordillera parent rocks such as slate (ca. 2,691 kg/m 3 ) and granite (ca. 2,750 kg/m 3 ) suggests that volume and density changes during weathering and erosion could have minor effects on volumetric comparisons between rock eroded from the fold-thrust belt and sediment deposited in the Subandean Zone.

In Eastern Cordillera, Interandean Zone, and Chaco Foreland
Using the eroded area from the balanced cross-section of Anderson et al. (2017) and the distribution of ages from various thermochronologic datasets used to infer timing of erosional exhumation (Ege et al., 2007;Anderson et al., 2018), the volume of material eroded from the Eastern Cordillera can be calculated for two time intervals. Across the Eastern Cordillera, rapid cooling was well underway by ca. 25 Ma and all apatite (U-Th)/He results show cooling below closure temperature by 15-12 Ma (Ege et al., 2007;Anderson et al., 2018). However, final exhumation of these samples may not have occurred until ca. 10 Ma development of the San Juan del Oro surface. Based on reasonable estimates of the paleogeothermal gradient in the region (e.g., Ege et al., 2007;Anderson et al., 2018), samples could have been at depths of up to ca. 2.5 km at 12 Ma, and subsequently brought to the surface by 10 Ma. Thus, between ca. 25 and 15-12 Ma, 133,818 km 3 was exhumed and removed from the Eastern Cordillera (Figure 6). By 10 Ma, an additional 61,307 km 3 had been eroded, totaling 195,125 km 3 of material removed from the Eastern Cordillera.
In the Interandean Zone, Triassic and younger rocks are no longer preserved Anderson et al., 2017). However, reset apatite fission track (U-Th)/He ages for Permian and older rocks require burial heating by a >4 km thick section of Mesozoic-Cenozoic rocks by late Oligocene-early Miocene time (Ege et al., 2007;Anderson et al., 2018). The eroded volumes from the western and eastern portions were calculated separately because they have distinct cooling ages determined from thermochronological data. In the western Interandean Zone, 34,488 km 3 material was removed between 25 and 18 Ma. In the eastern Interandean Zone, 31,446 km 3 was removed between 11 and 8 Ma.
In the Subandean Zone, remnants of the Cenozoic basin are preserved in the synclinal valleys between eroded anticlinal ranges. Erosion of the frontal fold-thrust belt contributed 54,132 km 3 material. To the east, the Chaco foreland has primarily been a region of sediment accumulation, with very  Uba et al., 2007). The Emborozú site has limited chronostratigraphic constraints, and no attempt was made to include the under-or overlying strata. At the Villamontes site, the impact of pure sandstone and pure shale lithologies on backstripping results was assessed, with minor differences in total tectonic subsidence between the lithologies, and recognizing that the Villamontes section is not a single lithology. For the Emborozú site, only sandstone lithology was used. Note that rapid sediment accumulation was ongoing at ca. 11 Ma, and similar rates of sediment accumulation are not observed at the Villamontes succession until ca. 8 Ma and later.
little shortening and exhumation. Calculations show a minimum volume of ca. 75,542 km 3 of preserved Cenozoic basin fill. This volume occupies the foreland basin between the Andean deformation front (near Villamontes) to the distal pinchout of the foreland basin projected ca. 190 km to the east (Figure 6).

From Beneath the San Juan del Oro Surface
After the preserved remnants of the San Juan del Oro surface were mapped ( Figure 7A) and interpolated to approximate the original geometry of the surface (Figure 7B), differencing the interpolated surface and the MERIT DEM yielded a ca. 96.5 km 3 FIGURE 6 | Gray polygons showing the areas that represent eroded material from different portions of the cross section. The areas were multiplied by the north-south extent of the map in Figure 1 to determine the eroded sediment volumes shown in this figure. Note that the Eastern Cordillera shows estimates for two distinct time intervals, and the Interandean Zone differentiates material eroded from the western and eastern segments. The basal dip of the foreland was continued eastward to determine the approximate position of the pinchout of the foredeep.
volume of material eroded from beneath the San Juan del Oro surface ( Figure 7C). The mapping of the San Juan del Oro surface is consistent with previous workers (Kennan et al., 1997;Anderson et al., 2018). These results show that the remnants of the surface today are not consistent with a horizontal or flat surface, but instead, an undulatory surface with a component of northward plunge. Whether this irregular surface represents the original geometry, or reflects subtle folding and deformation after the surface formed, remains unclear. Regardless, material removed from river valleys and gorges incised into the San Juan del Oro surface ( Figure 7C, black polygons) provides a minimum estimate of material supplied to the foreland since final construction of the surface.

Flexural Modeling
Five of six scenarios yielded viable flexural models (Figure 4). Only the scenario with linearly changing EET and forelandsloping topography did not yield any model results. In scenarios with constant EET (Figures 4A,B), modeled EET varied between 55 and 95 km. In the linearly changing case (Figure 4C), the left (west) side of the model had EET between ca. 10-50 km, and these increase to ca. 50-95 km on the right (east) side of the model. In the cases of exponentially changing EET (Figures 4D,E), the left (west) side of the model predicts a range of 5-90 km, and ca. 50-90 km on the right (west) side of the model.

Foreland Basin Subsidence
Results from magnetostratigraphic and backstripping analyses confirm that the Chaco foreland basin experienced an increase in sediment accumulation rate between ca. 11 and 8 Ma, recording ca. 1,750 m of tectonic subsidence. The parameters tested with the flexural models were assigned a large range of potential values in order to minimize bias in possible solutions. The purpose of these models was not to predict the exact values of specific parameters governing the evolution of the Chaco foreland basin, but to explore the range of conditions that could have generated the observed sediment accumulation patterns. The models are simplified (Figure 4), as they simulate emplacement of a load instantaneously, whereas the Eastern Cordillera and Interandean Zone developed over 10s of Myr, and have been translated hundreds of kilometers eastward during protracted shortening and underthrusting of the Brazilian craton. These results show that a range of load configurations could generate basin profiles outlined by backstripping and tectonic subsidence analyses. Likewise, the results show that a broad range of lithospheric conditions could have produced the observed magnitudes of tectonic subsidence (Figures 2, 5). The broad range of viable solutions gives support to previous interpretations of the Chaco foreland as a flexural foreland basin (e.g., Horton and DeCelles, 1997).
Given the dimensions of the Eastern Cordillera, some flexural model solutions are geologically unreasonable. For example, Figure 4A shows that one modeled solution (the narrowest, highest rectangular loads that could generate viable flexural profile) is unrealistically high (>7 km) and narrow (<100 km), because it is inconsistent with dimensions of the Eastern Cordillera.
It remains challenging to reconstruct ancient lateral changes in effective elastic thickness. However, these results highlight that left to right (west to east) increases in EET provide viable model solutions, consistent with the west to east increase in EET predicted across the western margin of South America (Tassara et al., 2007), and predicted for other segments of the Andes (Fosdick et al., 2014), although the majority of EET values predicted here are greater than those predicted by Tassara et al. (2007). The results from linear and exponential changes EET could be further explored to refine the paired shortening and subsidence histories of the Bolivian Andes. Regardless of remaining challenges with modeling and reconstructions, these results demonstrate that the tectonic subsidence observed in the southern Bolivian Chaco foreland are consistent with flexural subsidence predicted from a range of plausible Eastern Cordillera and Interandean Zone topographic loads.
The calculations of eroded volumes (Figure 7) in the study area support a genetic link between topographic loading and flexural subsidence. The amount of material removed from the Eastern Cordillera is 195,125 km 3 , which occurred before development of the high-elevation San Juan del Oro geomorphic surface. The sum of material eroded from the Interandean Zone and Subandean Zone is 119,766 km 3 . The volume of material removed from beneath the San Juan del Oro surface is 96.5 km 3 . The Chaco foreland contains at least 75,542 km 3 . These results are consistent with sediment provenance analyses showing that the fold-thrust belt was the primary source of sediment , and was capable of supplying all of the sediment in the Chaco foreland basin, even if the precise balance of material supplied westward to the Altiplano remains unconstrained. Barnes and Heins (2009) investigated a broader region spanning the Plio-Quaternary to recent, and demonstrated that incision into the San Juan del Oro surface provided 40-60% of Chaco foreland basin fill since ca. 2.1 Ma. These results show that there was ample additional sediment available from erosion of the fold-thrust belt to fill the Chaco basin since the middle Miocene.

Reconciling Foreland Basin and Fold-Thrust Belt Timing
Despite agreement between observed magnitudes of tectonic subsidence in the Chaco foreland (Figures 2, 5) and basin subsidence profiles predicted from various flexural models (Figure 4), and new recognition of ample detritus supplied from the thrust belt to the foreland basin (Figures 6, 7), there is an apparent discrepancy between the timing of exhumation in the fold-thrust belt beginning in the Eocene and the onset of rapid subsidence in the Subandean foreland basin in the Miocene (Figure 8). The incomplete Emborozú section shows rapid sediment accumulation was ongoing at 11 Ma , whereas the Villamontes section, situated eastward of the Emborozú section, experienced the shift to rapid sediment accumulation by ca. 8 Ma (Uba et al., 2007; this study) (Figure 2). This is consistent with eastward migration of a flexural foreland basin in front of an advancing fold-thrust belt.
We reassessed the thermochronologic results of Anderson et al. (2018) to compare exhumation in the fold-thrust belt with the known timing of rapid foreland basin sedimentation in the Chaco foreland. Broad trends in the thermochronologic data are consistent with the interpretations from Anderson et al. (2018) (Figure 8): after initial shortening along the axis of the Eastern Cordillera at ca. 40-35 Ma, which was focused just west of the Camargo syncline (Figures 1, 8), shortening propagated both westward and eastward. Westward (trenchward) verging contraction was accommodated by the Central Andean Backthrust Belt, which advanced to the Altiplano margin by ca. 27 Ma (San Vicente Thrust). Eastward (foreland) verging contraction advanced generally through the Cuesta de Sama anticlinorium, which may have been rapidly exhuming by 30-32 Ma, into the Interandean Zone by ca. 25-20 Ma, and into the Subandean Zone after 10 Ma. In this scenario, the limited chronostratigraphic results for foreland basin fill of the Camargo syncline region fall within the accepted range of possible depositional ages (between 55 and 25 Ma; DeCelles and Horton, 2003;Horton, 2005). Specifically, deposition near the Camargo syncline overlapped with early Eastern Cordillera shortening at ca. 40 Ma, and ceased by ca. 30 Ma when the Cuesta de Sama anticlinorium began to be uplifted east of the Camargo syncline depocenter. Deposition of the 2-3 km Camargo syncline succession over this time frame is consistent with the rapid accumulation rates observed in other proximal foreland basins.
In this reassessment of thermochronologic results, we highlight the assortment of apatite (U-Th)/He ages in the Eastern Cordillera. While higher temperature thermochronometers support westward advance of the backthrust belt, apatite (U-Th)/He ages from this same area have a range of reproducibility. Some samples show age uncertainties >15 Myr, while others show low uncertainties comparable to apatite (U-Th)/He results from the Interandean and Subandean Zones. This region of the Eastern Cordillera with higher age uncertainties is also where the San Juan del Oro surface developed. The contrast in style and magnitude of apatite (U-Th)/He thermochronologic cooling signals between the Eastern Cordillera and Interandean/Subandean Zones reflects a dichotomy of processes active during the Miocene across the foldthrust belt. In the Eastern Cordillera, the varied reproducibility among apatite (U-Th)/He samples occurs irrespective of mapped structures or proposed basement structural geometries. We suggest that the character of the apatite (U-Th)/He results in the Eastern Cordillera is more consistent with sedimentary processes, such as the infilling of intermontane basins and formation of geomorphic surfaces. Final construction of the San Juan del Oro surface was complete by ca. 10 Ma, but it may have begun developing earlier, after ca. 25 Ma cessation of most upper-crustal shortening in the Eastern Cordillera. Although the precise formation mechanisms of such high-elevation low-relief hinterland surfaces remain unclear (e.g., Fox et al., 2020) we envision a combination of erosional beveling, surface leveling, local deposition, relief reduction, and overall topographic homogenization. The combination of such processes would result in irregular cooling patterns expressed by thermochronologic datasets, attributable to alternating periods of localized erosion and deposition as topography was smoothed. Additionally, the structural restoration predicts lateral translation of this region over a long basement flat (Anderson et al., 2017), arguing against a clear structural control on cooling ages. The combination of a range of possible mechanisms associated with formation of the San Juan del Oro surface, and structural restorations, are consistent with the apatite (U-Th)/He results from the region. Thermokinematic modeling in central and northern Bolivia illustrate rapid cooling is focused above basement ramps, whereas cooling ages above basement flats are more variable because material is passively transported and exhumation rates are much lower (e.g., Rak et al., 2017;Buford Parks and McQuarrie, 2019).
In contrast with the Eastern Cordillera and San Juan del Oro region, apatite (U-Th)/He (and other thermochronometric) results from the Interandean and Subandean Zones display lower uncertainty and a clear eastward younging age signatures. This is consistent with eastward advance of Neogene shortening while the Eastern Cordillera was experiencing minimal upper-crustal shortening, developing the San Juan del Oro surface, and being transported along basement flats. The apatite (U-Th)/He results from the Eastern Cordillera suggest partial resetting during shallow burial but preclude the deep burial that might be expected from the volume of material known to have been removed from the Eastern Cordillera since 25 Ma. We are left with a questionwhere was the sediment eroded from the Eastern Cordillera and Interandean Zone being stored, if it was not sequestered in the Eastern Cordillera or deposited in the Subandean Zone before 11 Ma?
Addressing this question also requires satisfying other observations: (1) bivergent (westward and eastward) shortening in the Eastern Cordillera/Interandean zone fold-thrust belt that spanned ca. 40 and 10 Ma; (2) a contrast in apatite (U-Th)/He results from the Eastern Cordillera and Interandean/Subandean Zones suggestive of different processes acting in these regions to yield contrasting cooling ages; (3) combined investigations of eroded volumes, flexural modeling, magnetostratigraphic results, and backstripping that confirm that the fold-thrust belt provided sufficient sediment and topographic loading to explain observed Subandean basin volumes, despite (4) a temporal lag between the main phase of Eastern Cordillera/Interandean Zone shortening and subsidence in the adjacent foreland.
The following scenario is proposed. During initial shortening in the Eastern Cordillera, and until contraction reached the Cuesta de Sama anticlinorium by 32-30 Ma, structural restorations predict a >291 km distance between the deformation front and the Villamontes section. The Eocene foreland basin that developed in the Eastern Cordillera and Interandean Zone (best preserved in the Camargo syncline; DeCelles and Horton, 2003) likely received a significant volume of the material eroded from the axial and western segments of the Eastern Cordillera. During this time, these segments of the Eastern Cordillera were eroded to Paleozoic levels, such that subsequent Cenozoic intermontane basin fill was deposited on Ordovician strata in angular unconformity (Horton, 1998(Horton, , 2005, and flanked by development of the San Juan del Oro surface between ca. 25 and 10 Ma. While the initial foreland basin was accumulating sediment, the future Subandean zone was not. The Subandes may not have received sediment at this time even in the backbulge depozone if the flexural wave was damped out (DeCelles and Horton, 2003). Passage of the forebulge depozone, recorded as paleosols in the Oligocene to middle Miocene Petaca formation, suggest that the Subandes were not in the foredeep depozone in the Eocene. To constrain the basin width and lithospheric parameters, we assume that the distance to the Villamontes section, 291 km, represents the maximum foredeep width. Assuming Young's modulus of 70 GPa, Poisson's ratio of 0.25, mantle and crust density of 3,300 and 2,600 kg/m 3 , respectively, an effective elastic thickness >40 km would have been required to generate a foredeep depozone ca. 291 km wide. Tassara et al. (2007) report the effective elastic thickness of the modern lithosphere under the Andean orogen to be generally less than 15 km. This suggests that the effective elastic thickness of the lithosphere during initial shortening would have been between 40 and 15 km. Even as shortening advanced eastward, the Villamontes locality remained far to the east, beyond the flexural wavelength of the foredeep depozone basin, until the late Miocene. The Villamontes succession would have entered the rapidly subsiding foredeep depozone only after significant advance of the foldthrust belt and underthrusting of the Brazilian Craton, which would have decreased the distance between the Eastern Cordillera topographic load and the Subandean zone, and increased the effective elastic strength of the lithosphere, making the depozone wider. Between 40 and 10 Ma, material eroded from the foldthrust belt was continuously recycled in the advancing foreland basin, as demonstrated by sediment provenance results . Since 10 Ma, additional erosion of the Subandean Zone, and incision into the San Juan del Oro surface further supplemented sediment delivery to the basin. A migrating foreland basin system is consistent with past interpretations of an advancing flexural load and progradational depositional systems (e.g., Uba et al., 2006). This work is consistent with previous findings (DeCelles and Horton, 2003;Uba et al., 2006), and reconciles apparent temporal inconsistencies among the timing and position of shortening, eroded sediment volume and sediment accumulation timing, and geomorphic observations.
In addition to the clear tectonic mechanisms the controlled foreland basin evolution in the central Andes, Miocene climate shifts must also be considered. Various changes in depositional setting and style at ca. 8 Ma, including a switch from braided streams to a progradational fluvial megafan setting in the Villamontes section, led Uba et al. (2005) to postulate that climate change played a dominant role in modifying foreland sedimentation rates. Strecker et al. (2007) stated that the South American monsoon began between 12 and 8 Ma, and that the central Andean region changed from a semi-arid or arid system to a humid environment by middle to late Miocene time. The links among shortening, erosion, and foreland basin sedimentation discussed in this manuscript provide a mechanism to explain the shifts in sedimentation rate apart from climate. A forelandyounging pattern of rapid sediment accumulation is revealed by comparison of separate localities (the Camargo syncline, Emborozú, and Villamontes sections). This time-transgressive history is consistent with a continuous eastward advance of the fold-thrust belt, and inconsistent with a single shift in sediment accumulation rate by a climate-driven process. We envision that Miocene climate shifts modified sediment erosion, transport, and depositional processes that further amplified tectonic processes. These include protracted fold-thrust belt advance toward the foreland, exhumation and sediment recycling, and flexural subsidence, which acted on development of the Chaco foreland basin. The results presented here bolster the recognition that interactions between climate, tectonic, and surface processes on orogen development are key, despite persistent challenges disentangling their precise contributions.

CONCLUSION
This study investigated the linkages among shortening, erosion, sediment transport, and foreland basin subsidence in the southern Bolivian Andes. Numerous datasets were brought to bear on the problem: structural geometries and exhumation timing provided by thermochronologic datasets spanning the Altiplano to Subandean Zone; new magnetostratigraphic results from a Subandean section at the Villamontes succession to generate a detailed chronostratigraphic framework; backstripping of the Villamontes and Emborozú sections to determine the magnitude of tectonic subsidence and quantify rates of sediment accumulation; flexural modeling to investigate the load dimensions and lithospheric parameters that accommodated observed foreland basin geometries; and eroded volume calculations to determine the magnitude of material removed from different sections of the Andean fold-thrust belt, and the minimum volume of sediment stored in the modern Chaco foreland.
This work shows that the magnitude of foreland basin sediment accumulation is consistent with reasonable fold-thrust belt geometries and lithospheric conditions. Further, the foldthrust belt would have been the dominant sediment source by volume and provided sufficient material to fill the foreland basin. Material excavated from beneath the San Juan del Oro surface contributed minor sediment volume to the foreland basin. Despite these consistencies, a key discrepancy remains. Specifically, how the modern foreland basin could display middle Miocene and younger onset of rapid sediment accumulation, despite evidence for fold-thrust belt shortening and erosion since ca. 40 Ma. We propose that the large magnitude of shortening and attendant lateral translation of the fold-thrust belt, coupled with flexural subsidence over a weak lithosphere, created an initial foreland basin that was too narrow to induce sediment accumulation in the distal Subandean Zone. The early flexural foreland basin has been largely erosionally recycled as the foldthrust belt advanced eastward. Remnants of this early foreland basin have been identified in the eastern segment of the Eastern Cordillera (Camargo syncline and related sections). Only in the past ca. 11 Myr has a combination of crustal shortening, attendant underthrusting of the Brazilian Shield, and foreland advance of the fold-thrust belt been sufficient to bring the Subandean localities (Villamontes and other sections) into the proximal foredeep depozone. An interesting possibility emerges, such that once the early foreland basin is entirely eroded away, this segment of the Andes may have a foreland basin stratigraphic record that does not preserve the initial 10s of Myrs of early Andean orogenesis.

DATA AVAILABILITY STATEMENT
The original contributions generated for this study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
All authors were involved in the conception of this manuscript. NP and BH conducted the magnetostratigraphic analyses. NP conducted the backstripping, flexural modeling, most figure creation, and most manuscript writing. RA calculated the eroded areas from the cross section. BO conducted the San Juan del Oro eroded volume work.

ACKNOWLEDGMENTS
Magnetostratigraphic data in this study comes from the undergraduate thesis of NP, who thanks undergraduate honors program mentors Chris Bell, Bill Carlson, Jack Holt, Isaac Smith, and John Kappelman from the UT Austin Paleomagnetics Laboratory. NP and BH especially thank Cornelius Uba for leading us through Cenozoic sections in Bolivia, field discussions, and assistance with magnetostratigraphic sampling. BO thanks George Allen for assistance with ArcGIS methods.