Structural Style and Kinematic History of the Colombian Eastern Cordillera

Fold-and-thrust belts and their associated structures are among the most common geological features of convergent margins. They provide significant information about crustal shortening and mountain-building processes. In subaerial belts, where the erosional rates are high and the growth strata are mostly eroded, methodologies such as that presented here can provide insights into to their formation. Two 2D cross-sections located in the Eastern Cordillera of Colombia are presented in this research. These sections extend from the Bogota Savanna to The Llanos, parallel to the regional deformation direction. Section construction was carried out using commercial surface data, and seismic information provided by Ecopetrol. Published thermochronometric data, gravel-clast petrography analysis, and paleoflora analysis were used to construct a viable tectono-evolutionary history of the study area. This evolutionary model is presented here in two palinpastic restorations from the Early Paleogene to Recent (∼65 Ma to Present-day). Section 1 and Section 10 accumulated 17.3 km and 19.5 km of shortening, respectively. The section reconstruction displays two major tectonic events – post-rift subsidence during the Early-Mid Paleogene, and positive inversion from the Oligocene to Recent (∼33 Ma to Present-day). This investigation focuses on the compressional period, where the structural analysis evidences an acceleration in the shortening rate, as well as a progressive migration of the deformation from northwest to southeast. This research discusses the extent and limitation of this methodology, as well as the principal structural aspects of the reconstruction.


INTRODUCTION
Fold-and-thrust belts are the typical place in which shortening is accommodated in the crust, and they are widely distributed (Nemcok et al., 2005;Cooper, 2007). For hydrocarbon exploration, the knowledge of the structural evolution of the fold-and-thrust bels is critical to understand the generation, migration, and accumulation of oil and gas (e.g., Masini et al., 2011).
Two balanced cross-sections were constructed in the Eastern Cordillera of Colombia based on detailed maps, surface and well dips, and confidential seismic data. These cross-sections are in the east portion of the Eastern Cordillera, between latitude 5.60°N and 4.00°N, and longitude 74.30°W and 72.60°W (Figure 1). For dating thrust events, a wide range of methods can be applied. The timing of thrusting, and associated folding, is best determined when pre-kinematic, syn-kinematic, and post-kinematic depositional sequences are preserved in the study area. Furthermore, when the syn-kinematic deposits (i.e., growth strata) are well recorded, not only the timing, but also the folding mechanism can be recognized in detail (Suppe and Medwedeff, 1990;Erslev, 1991;Poblet et al., 1997;Shaw et al., 2005). This approach has been used in numerous studies in deep-water fold-and-thrust belts (e.g., Bilotti et al., 2005;Rowan and Peel, 2005;Vidal-Royo et al., 2013;Butler, 2020). However, in subaerial belts, where the erosional rates are high, and the growth strata are mostly eroded, as in the Eastern Cordillera, different methods can provide insight into the structural evolution.
Both balanced cross-sections were restored from the Early Paleogene to Recent (∼65 Ma to Present-day). The kinematic evolution of these cross-sections was time-calibrated by means of thermochronometric data, and the chronostratigraphy of the syntectonic successions shaped by the contractional tectonic regime. The necessary topographic surfaces for each restoration stages were inferred from published paleoflora analysis (Wijninga, 1996;Hooghiemstra et al., 2006;Anderson et al., 2016). The structural analysis of these two cross-sections allowed for the understanding of the strain distribution and shortening rate in the study area.
This study focuses on the structural evolution of a portion of the Eastern Cordillera by a comprehensive and multidisciplinary methodology. It analyses the extent and limitation of this methodology, compares the most relevant results with previous researches, and discusses the principal structural aspects of the compressional period in this area (from the Oligocene to Recent).

GEOLOGICAL SETTING
Structural Setting the inversion of an earlier rift system (involving basement uplift), and a subsequent activation of footwall shortcuts. Inversion onset could be established in the Late Cretaceous , with an increase in the compressive activity from the Pliocene to Recent, and the active deformation continues to the Present-day .
The study area ( Figure 1B) is located between latitude 5.60°N and 4.00°N, and longitude 74.30°W and 72.60°W (Eastern flank of the Eastern Cordillera). Along the dip direction, this part of the Cordillera displays overall symmetrical deformation (Figure 2), and can be divided into three regions, the Western and Eastern marginal portions, and the Bogota Savanna in the center. The Western and Eastern margins appear to have accommodated the shortening in a brittle manner. In these areas, the inverted structures involve basement uplift that gave rise to thickskinned systems ( Figure 2). In contrast, the most marginal areas exhibit thin-skinned tectonics and they are bounded to the East by The Llanos foreland basin and Middle Magdalena Valley intra-mountain basin to the West ( Figure 1A and Figure 2). The Bogota Savanna (Figures 1, 2) is located in the highland between the marginal fold-and-thrust belts and is characterized by a low relief that is approximately 2600 m above sea level . Unlike the marginal systems, the deformation along this area features shorterwavelength folding with low amplitude. The possible causes of this distinctive deformation are discussed below.

Stratigraphy
The stratigraphy in the area (Figure 3) can be broadly divided into Paleozoic, Mesozoic, and Cenozoic rocks depicting multiple regional unconformities (Martinez, 2006).
Describing the stratigraphic column from older to younger (Figure 3), the oldest drilled Paleozoic rock is Ordovician (Martinez, 2006). In this study, the Paleozoic rocks are considered the economic basin basement, and the exhumation of this unit resulted in up to 4 km of low and medium grade metamorphic rocks and intermediate to acid intrusives (Parra et al., 2009b;Segovia, 1965). The Devonian to Carboniferous rocks are deposited disconformably overlying the basement, and they can reach up to 4 km of thickness. They are formed by clastic platformal sequences described as the Devonian-Carboniferous Farallones Group ( Figure 2) located in the Quetame Massif (Ulloa and Rodriguez, 1976).
During the early period of rifting, from the Triassic to Jurassic, the deposition was predominantly volcanoclastic with some intercalation of non-marine to shallow marine sediments (Parra et al., 2009b). These sediments were locally deposited on the narrow asymmetric grabens located mostly in the western flank of the Eastern Cordillera (e.g., Rusia and Giron Formations in Figure 3).
The thickness (approximately 4-6 km) and wide distribution of the Lower Cretaceous rocks, suggest that the Later Mesozoic rifting phase initiated a wider extensional system (Mora et al., 2006;Mora et al., 2009). During this period the basin was bounded by regional extensional faults dipping basinward (La Salina fault to the west and Servita-Lengupa system to the east in Figure 1B). Macanal, Las Juntas, and Fomeque Formations were deposited in the study area during this period ( Figure 3).
During the post-rift and thermal subsidence (Sarmiento-Rojas et al., 2006), in the Upper Cretaceous, approximately 1.5-2.0 km thick sediments were deposited in this basin. They are Une, Chipaque, and Guaduas Formations and Guadalupe Group ( Figure 3).
The Late Cretaceous inversion of this Mesozoic rift basin gave rise to a foreland system east of the basin which extends throughout the foothill and The Llanos basin (Bayona et al., 2008;Cooper et al., 1995). The Cenozoic rocks in the Medina Basin ( Figure 1B) comprise a Paleocene to Late Miocene sequence which evolves from coastal plain and tidally influence lacustrine sediments of the Barco, Los Cuervos, Mirador, Carbonera, and Leon Formations (Figure 3), to the proximal and alluvial deposits of Guayabo Group (Cooper et al., 1995;Parra et al., 2009b).

METHODOLOGY AND DATA
The methodology applied in this study covers a number of structural-geology subjects, such as section construction and FIGURE 2 | Regional 2D structural cross-section along the Eastern Cordillera of Colombia (modified from Teixell et al., 2015). It is characterized by an overall symmetrical deformation with two marginal systems and a highland in the middle of the section (Bogota Savanna). The inversion-tectonic portion forms a thick-skinned system, whereas the foothill and foreland basin have experienced a thin-skinned deformation.
Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 636458 balancing, low-temperature thermochronology-based section restoration, and paleotopograhy reconstruction. For this study, Ecopetrol Colombia has provided commercial datasets: 2D seismic reflection, surface dips, geological maps, and well data (formation markers and strata dip). The thermochronology used to constrain the structural evolution has been published by Parra et al. (2009a). The paleotopography construction was carried out based on the paleoaltimetry and paleoflora analysis published by Van der Hammen et al. (1973), Hooghiemstra et al. (2006), and Wijninga (1996); and gravel-clast petrography analysis from the Medina Basin (Parra et al., 2009b;Parra et al., 2010). This comprehensive method can give insights into the evolution of the fold-and-thrust belts where the growth strata are highly eroded. However, the accuracy of the outcomes is subject to the quantity and quality of the available data.
study. The 'kink method' (cf. Allmendinger, 2015) was applied for the section construction and balancing. For the two sections presented in this study, the section construction was carried out by a comprehensive methodology which involves different data from diverse sources. Figure 4 displays an example of section construction and balancing, and the multi-source data involved in this process. The topographic profile on each section (Figure 4) was extracted from the Digital Elevation Model (DEM) available for this research (International Hydrographic Organization, 2003). The intersections of the different formations ( Figure 1B) with the topographic profile were calculated projecting the formational contacts of the geological map as displayed in Figure 4. Similar process was used for the fault intersection. Surface dips ( Figure 4) were used to constrain the strata dip angles nearby the topographic surface. Limited and confidential 2D seismic reflection and well data were used to interpret the deep structures of the foreland basin; particularly in the Medina Basin and The Llanos where thick layers of recent sediments preclude the use of surface-observation method.

Thermochronometric Dating
Techniques involving low-temperature thermochronometric dating have been widely used to reveal the exhumation, and its associated cooling history, in many geological settings affected by upper-crust processes (e.g., Garver et al., 1999;Reiners and Brandon, 2006;Blythe et al., 2007). This method allows a correlation between the rock exhumation and the processes that transport rocks to the topographic surface. Methods such as fission track (FT) and (U-Th) He have proven successful in areas where the exhumation occurred during the Late Mesozoic to Cenozoic period; encompassing extensional tectonic settings (e.g. Armstrong et al., 2003;Ehlers et al., 2003;Blythe et al., 2007), and compressional systems such as subduction zones and collisional belts (e.g., Garver and Kamp, 2002;Ehlers et al., 2003;Garver et al., 2005;Cao et al., 2013). Specifically, in the Easter Cordillera, this technique has allowed for significant advancement of understanding about how the Eastern Cordillera formed (e.g., Parra et al., 2009b;Mora et al., 2010;Parra et al., 2012;Mora A. S. et al., 2013;Cao et al., 2013;Carrillo et al., 2016).
Thermochronometric dating does not provide structuralevolution information per se, but the thermal history of an exhumed sample. However, the temperature variation of a sample that has been buried relates to the exhumation of this sample (Fitzgerald et al., 1995;Lisker et al., 2009). Moreover, exhumation can be described as the displacement of rocks with respect to the topographic surface . Therefore, exhumation cannot be simply associated with an uplift rate; but to a more sophisticated process -denudation. Erosion or tectonics (and a combination of both) generate denudation at the Earth's surface. Erosion is the process by which soil or rock experience breaking down, chemical solution and transport (Michael, 2008); whereas tectonic denudation is associated with extension and normal faulting resulting in a rapid removal of large rock volumes (Lisker et al., 2009). Due to the variety of processes associated with the thermal history of a sample, it is necessary to constrain this information with geological data to obtain a verisimilar solution in terms of structural evolution (i.e., uplift, folding, and transport rate). In this research the thermochronometric data from apatite and zircon fission tracks are combined with paleoelevation and paleoprovenance studies to sequentially restore the deformation in the study area.
The fission track method is based on the fission decay of uranium which forms a linear damage trace in the mineral lattice (Fleischer et al., 1975). The fission track traces form continuously at relative low temperature, but they are annealed rapidly if the mineral experiences high temperature. During the annealing process, the damage traces reduce the length until they vanish away; therefore, the previous thermal history is reset. Subsequently, if the annealed sample is brought to temperatures below the annealing temperature, fission tracks forms over again recording a new thermal history. For apatite, the temperature of continuous annealing (or closure temperature) is approximately 100-120°C, and for zircon it is approximately 250-300°C (Green et al., 1986). A relevant aspect to consider about this method is that there is a range of temperature of partial annealing, the "partial annealing zone" (PAZ). Within the PAZ the annealing process intensifies rapidly when the temperature increases. Conventional ranges of PAZ from ∼60-120°C and ∼210-300°C are estimated for apatite and zircon respectively (Green et al., 1986;Tagami, 2005). If a sample experiences slow transit through the PAZ, the annealing process would affect the tracks length precluding a confidence analysis. In orogenic belts, which have experienced rapid cooling due to uplift and erosion, such as the Eastern Cordillera, a characteristic fission-track profile should be recorded.
The geographic distribution of the apatite and zircon fissiontrack samples used for this study is shown in Figure 5. These 49 samples have been analyzed by Parra et al. (2009a) to study the diachronous exhumation of the Central and Eastern Cordillera. Table 1 presents the fission track ages with their corresponding formation. Further details regarding the sample analyses are presented in Parra et al. (2009a) Tables 1, 2, and its respective auxiliary material.

Paleotopographic Restoration
In this study the sequential structural restoration has been carried out by the analysis of the exhumation history of a series of fission-track samples. As mentioned above, exhumation can be defined as the displacement of a portion of rock with respect to the topographic surface . Thereby, the reconstruction of paleotopographies becomes crucial for the process of relating exhumation and cooling FIGURE 5 | Geological map with the apatite and zircon fission track sample distribution, and the trace of the two cross-sections restored in this study.
Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 636458 history to uplift and erosion. A paleotopography can be modeled by means of indirect methods, and can be subject to significant errors and uncertainty depending on the quantity and quality of the data involved in the process. In this research the paleotopographies were inferred based on published analysis of paleoaltimetry and paleoflora for the axial zone and Bogota Savanna (Wijninga, 1996;Hooghiemstra et al., 2006), and lipid biomarker analysis (Anderson et al., 2015;Anderson et al., 2016). Paleoaltimetry and paleoflora: these methods have been used previously in the Eastern Cordillera to gain understanding about the structural evolution of this orogen (e.g., Hooghiemstra and Cleef, 1995;Hooghiemstra et al., 2006;Torres et al., 2005;Veer and Hooghiemstra, 2000). This  methodology allows the correlation between fossil pollen spectra with similar pollen spectra from modern vegetation. Figure 6 summarizes the fossil flora analysis conducted by Hooghiemstra et al. (2006), and Wijninga (1996), in the Late Miocene to Late Pliocene fluvialite Tilata Formation, from samples collected in the Bogota Savanna. They argue that for the Middle Miocene the paleoelevation in the Bogota Savanna region was ∼700 ± 400 m, during the early Pliocene ∼1100 ± 400 m, in the Late Pliocene ∼2000 ± 400 m, and for the Quaternary ∼2500 m -which is the Present-day average of elevation in this region. This information has been previously used in the same fashion in the Eastern Cordillera to calibrate thermochronometric restorations (e.g., Wijninga, 1996;Mora A. S. et al., 2013;Carrillo et al., 2016). More recently, the chronology of deposition of those units and the paleoelevaton estimations have been reassessed by Anderson et al. (2015); however, the results essentially confirm those by Wijninga (1996).

Provenance Analysis and Unroofing History
The gravel-clast petrography analysis allows the temporal and spatial correlation of eroded and transported sediments with their final depositional area. Foreland basins are generally formed as a result of flexural subsidence which is driven by the thrust-sheet loading affecting the frontal system (DeCelles and Giles, 1996). Therefore, the sediments deposited in the accommodation space of the foreland basin derive principally from the frontal-system erosion (Dickinson and Suczek, 1979;DeCelles and Giles, 1996). Hence, the gravel-clast petrography analysis of a foreland basin permits the reconstruction of the erosion windows on the paleoprovenance areas. Parra et al. (2010) tracked the occurrence of different conglomerate clasts deposited in the Medina Basin during the Oligocene-Miocene to evaluate the unroofing history of the exposed Quetame massif (Figure 7). They interpreted a normal unroofing sequence suggesting that Paleogene rocks are the sources of the conglomerates below the Upper C1 Member; and the sources of the Upper C1 Member, Leon, and Guayabo Formations are Upper Cretaceous rocks (Guadalupe, Chipaque, and Une Formations). In this research, this information is used to calibrate the paleotopography reconstruction in the frontal system and adjacent areas (Quetame massif, Medina Basin, and the proximal part of The Llanos). Even more, it is necessary to recreate the erosion windows in the frontal system that contributed to the depositional record in the Medina Basin.

Paleocurrents and the Age of the Guavio Anticline
The Guavio anticline is a thin-skin structure that is bounded by the Guaicaramo fault to the southeast ( Figure 4). Seismic data suggest that the Guavio anticline is a fault bend fold related to the subsurface shape of Guaicaramo fault. If that is the case, the age of the Guavio anticline would be directly associated with the age of the Guaicaramo fault. This fault cuts through the Neogene and Quaternary sedimentary units in The Llanos foreland ( Figure 4). There are also active terraces related to this fault . Therefore, it should have been active most likely from the latest Neogene throughout the Quaternary. However, Quintero (2010) documented a significant shift in paleocurrents in the Neogene sedimentary units in the Guavio anticline, from dominantly parallel to the structural grain to perpendicular FIGURE 6 | Paleoaltimetry estimations based on fossil pollen spectra analysis for the Bogota Savanna and surrounding areas (modified from Hooghiemstra et al., 2006;Wijninga, 1996).
Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 636458 and pointing to the NW. This shift was observed by Quintero (2010) in the uppermost Miocene units. Therefore, the most evident growth of the Guavio anticline can be bracketed between Late Miocene and the Present-day. Uncertainties would be related to the age of the sedimentary units where the paleocurrents were measured (Quintero, 2010).

2D Thermokinematic Restoration
The integration of the above-mentioned techniques and data sets was used to constrain the restoration at different times for both cross-sections. The restorations have to respect different assumptions like the following. Thermochronology would guide at which temperature the different rock units should be for the restoration steps. While this happens, the unroofing history from gravel petrography would limit which units are exposed and shedding detritus at the surface. Therefore, the still buried and hot rocks in the subsurface (and their thicknesses) should be coherent with the actively eroding and exposed rock units and, more important, should not be the same. For example, while the actively eroding units at a particular time are at the surface, thermochronometers may predict that other units should be 4 km underneath the topography. Thus, the calculated thicknesses for the different rock units, which is still buried and preserved, should respect both data sets.
Finally, the paleoelevation estimates would allow to approximately constrain the maximum elevation above the sea level where the exposed rock units could be. This should be in agreement with the predicted fault geometries at depth, and the structural and topographic elevation where the rocks should be exposed. Thus, the above-mentioned data relate denudation, rocks, and surface uplift in a single structural restoration. If the available data for this study area are carefully used, the deformation of the different rock units would be constrained with reasonable accuracy.

RESULTS
The 2D section construction, balancing (Figures 8,9), and restoration ( Figure 10) of Section 1 and 10 are discussed here. It was possible to obtain nine evolutionary stages from the thermochronology constraints. These evolutionary stages can be broadly divided in two major tectonic events: post-rift, and inversion tectonics. The Cretaceous rift period was followed by a post-rift regime characterized by cooling and subsidence. From ∼65 Ma to ∼40 Ma (Early Paleocene to Late Eocene) the post-rift tectonics facilitated the accommodation space for the deposition of the following Paleogene units: Guaduas Formation at 65 Ma ( Figures 10A,B), Los Cuevos Formation at ∼55 Ma ( Figures 10C,D), Mirador Formation at

Structural Style
One of the principal structural features of this area is the coexistence of thick and thin-skinned deformations (Figures 8,9). In the Eastern Cordillera, from the Servita-Lengupa fault toward the hinterland, a thick-skinned deformation with basement uplift caused by the inversion tectonic initiated during the Oligocene is observed, and also documented in previous works (e.g., Horton et al., 2010;Ramirez-Arias et al., 2012). Farther to the Southeast, in Medina Basin (Figure 9), the occurrence of shallower detachments associated with thrust ramps results in a thin-skinned deformation of more recent age. The hinterland of the study area is the locus of a high plateau (cf. Ramirez-Arias et al., 2012), known as Bogota Savanna (Figures 8,9). Section 1 and 10 run along the middle and northeast of this structural elevation respectively ( Figure 5). Paradoxically, the degree of deformation on this plateau is significantly less than in the foreland system; however, the plateau is in a much higher structural position. A possible cause of the plateau uplift is discussed below.
From the Bogota Savanna toward the fontal system ( Figures  8, 9) the most important structural and topographic relief is observed. This deformation is associated with the external inversion faults (Servita-Lengupa, and Tesalia faults) which accommodate the thick-skinned displacement. Farther to the foreland (i.e., to the southeast) the foreland basin is disrupted by a series of low-angle (∼30°) thrust faults (Mirador and Boa fault in Section 1, and Guaicaramo fault in Section 10) which give rise to a thin-skinned system. They are interpreted as footwall edge shortcuts (cf. Rowan and Linares, 2000;Tesón et al., 2013) connected to the aforementioned inversion faults (Figures 8, 9). In the case of the Guaicaramo fault, it can be described as a multi-ramp thrust with a sub-horizontal segment in between. The Medina (or Guavio) anticline in the foreland basin in Section 10 ( Figure 9) is the consequence of a fault-bend-folding deformation  The deformation process of the inversion tectonics from the Early Oligocene to Recent on both section is described in the following.

Configuration Prior to Shortening ca. 33 Ma (Early Oligocene)
Previous studies based on thermochronological information, stratigraphic sections, and sedimentary petrology (e.g., Horton et al., 2010;Parra et al., 2009b;Parra et al., 2010;Ramirez-Arias et al., 2012) present evidence which date the initial exhumation in the Eastern Cordillera during the Oligocene. Hence, the ∼33 Ma restoration (Figure 11) is considered the youngest structural setting before the compressional tectonics; also, the maximum sedimentary load is assumed for this stage. The C8 Member (of the Carbonera Formation) in the foothill and The Llanos, and the Concentracion Formation on the Eastern Cordillera were deposited during this age (Cooper et al., 1995). The C8 Member is described as marine-influenced and lower coastalplain sedimentation (Cooper et al., 1995), therefore the paleotopography is assumed to be a flat-lying surface vertically located near the mean sea level (Figure 11).
Assuming a thermal gradient of approximately 25°-30°C per kilometre in depth (Allen and Allen, 2006), the closure temperature of apatite fission tracks (∼120°C) can be estimated with a line which is parallel to but 4 km underneath the paleotopography. Applying the same reasoning, the closure temperature of the zircon fission track (∼240°C) can be estimated at 8 km beneath the paleotopography. This apatiteclosure-temperature approximation in all models is demarked by a light-blue dashed line and referred to as ACTA; whereas the zircon-closure-temperature approximation is referred as ZCTA, and is displayed as the light-red dashed line (Figures 11-16).
For the necessity of projecting the samples onto the sections, some thermochronometers display an offset with respect to the Present-day topography (in dashed red line in Figure 11). Additionally, some small amendments (vertical movements) were made on the projected samples to honor their relative position in relation to the formations from which they were extracted. To clarify the analysis of the thermochronological data both sections have been divided into three zones. Zone-1 covers the Foothill and the frontal thrust system of the Eastern Cordillera ( Figure 11). Zone-2 extends across the central portion of the sections, which is featured in the Present-day deformation by a marked monocline dipping northwest. Zone-3 covers the area of the Bogota Savanna. Note that Section 1 does not have thermochronometers located over zone-3, therefore the deformation and analysis of this zone are inferred from Section 10.
In zone-1 of Section 10, ZFTs Z05 and Z06, which are located in Macanal Formation ( Figure 11B) display an average age of 16 Ma. They are considered reset, and this is consistent with their relative location with respect to the ZCTA as they are displayed beneath the light red-dashed line. A similar analysis can be applied to sample Z07 which is located in the Farallones Group. From the AFT samples A27 and A28, it is also inferred that the samples are reset with and average age of 2.0 Ma, and located beneath the ACTA. In zone-1 of Section 1 ( Figure 11A) there are three ZFT samples, Z16 (in Quetame Group), Z19 and Z20 (both in Buenavista Fm). Z19 and Z20 display anomalous ages of 145.2 Ma and 165.9 Ma respectively. These dates markedly contrast with the overall knowledge about the structural evolution of this belt, therefore it is assumed that they have never reached the zircon closure temperature (∼240°). Consequently, they can be observed in Figure 11B above the ZCTA. In this zone the apatite fission track age is in an average of 2.7 Ma with all the samples being reset (A41, A44, A46, A45). Zone-2 of Section 10 ( Figure 11B) is populated by eleven AFT samples. Every sample appeared to be reset, apart from A13 which displays an anomalous apatite age of 38 Ma, indicating that it could be partially annealed. To the southeast portion of this zone samples A16 to A18, A25, and 26 are located in the Early Cretaceous formations and they have an average apatite age of 6.7 Ma. To the northwest part of this zone, samples A12, A14, A15, A22, and A23 are located in the post-rift (Late-Cretaceous) formations with an average apatite age of 8.4 Ma. In zone-2 Section 1 ( Figure 11A), the fission track samples were collected from the deepest stratigraphic units exposed on the Present-day topographic surface (from the pre Devonian Quetame Group to the Early Cretaceous formations). In the southeast portion of this zone the ZFT samples Z8 to Z10, Z15, Z17 and Z18 appear to be reset, with an average zircon age of 10.0 Ma. The average apatite age calculated from A29 to A34, A40, A42, and A43 in this area is 2.1 Ma. In the northwest portion of this zone two zircon fission track samples Z11 (136.5 Ma), and Z12 (61.3 Ma) appear not to be reset, therefore in the pre-deformation model they are located above the ZCTA. The apatite fission track samples A35 to A37, placed in the Lower Cretaceous formations to the northwest of this zone, have an average apatite age of 10.0 Ma.
In zone-3 of Section 10 ( Figure 11B) there are three AFT samples, A19 to A21. A19 is from the Paleogene Upper Socha Formation, having an apatite age of 45 Ma, which indicates that it could be partially annealed. A20 (17.6 Ma) and A21 (24.4 Ma), which are consider reset, are located in the Paleocene Cacho Formation and Late Cretaceous Guadalupe Group respectively.

Restoration ca. 23 Ma
At the beginning of the Miocene age (∼23 Ma) the C5 Member (of the Carbonera Formation) was deposited on the foreland basin ( Figure 12). The inferred paleoelevation analysis ( Figure 6) indicates that before the Middle Miocene (∼15 Ma) the paleotopography was similar to the mean sea level in the region. This is consistent with the marine influenced origin of C5 Member (Cooper et al., 1995).
The deformation initiated during the Middle Eocene to Oligocene (Parra et al., 2009b;Horton et al., 2010;Parra et al., 2010;Ramirez-Arias et al., 2012) along the axial zone. After the initiation, the western portion strain has propagated toward the Magdalena Valley, whereas the eastern portion toward The Llanos in an overall symmetric fashion ( Figure 2). The eastward deformation is reflected in this model with the initial deformation localized in zone-3 ( Figure 12). It is in zone-3 of Section 10 that the AFT samples with oldest apatite ages are extracted. The reset samples in this zone (A20 and A21) have an average apatite age of 20 Ma, therefore during this period they should be located near the ACTA -∼4 km in depth ( Figure 12B). The deformation in this zone, but in Section 1 ( Figure 12A), has been inferred by taking into account the lateral continuity from Section 10 toward the south. The accumulative shortenings from ∼33 Ma to ∼23 Ma, are 0.6 and 3.0 km for Section 1 and 10, respectively.
It has been mentioned above, that the onset of deformation estimated in the literature for the western structural domains (Saylor et al., 2012a;Saylor et al., 2012b), which is the zone-3, occurred by the Late Eocene. Therefore, the deformed section at 23 Ma actually represents the finite shortening that occurred between the Late Eocene and Late Oligocene. In this context, the available apatite ages give a minimum age for the exhumation and associated shortening.

Restoration to ca. 11 Ma
At the Late Miocene (∼11 Ma) the Leon Formation was deposited in the study area ( Figure 13). An interpolation of the paleoelevation information ( Figure 6) has been used to model the paleotopography. At the hinterland the elevation is assumed to be approximately 1000 m, progressively decreasing toward The Llanos until reaching similar values to the mean sea level.
During the period between ∼23 Ma and 11 Ma, in zone-3 on both sections, the exhumation continued with the AFT samples in a cooler position above the ACTA. The deformation propagated toward zone-2 and northwest of zone-1. The structural evolution from ∼23 Ma to ∼11 Ma is featured in the inversion of the Servita-Lengupa fault system forming a gentle and long-wave anticline ( Figure 13). This anticline is asymmetric in Section 10 with a lengthy back limb dipping toward the northwest; whereas in Section 1 it appears to be symmetric ( Figure 13).
Every ZFT sample located along Section 10 was in a cooler position with respect to the ZCTA ( Figure 13B). In Section 1, at the southeast part of zone-2, ZFT samples (Z08 to Z10, Z15 to Z18) with an average zircon age of 9.9 Ma were distributed near the ZCTA ( Figure 13A). This constrains the anticline symmetry since it generates a steeper and shorter backlimb compared to the anticline geometry observed on Section 10. To the northwest part of zone-2 the reset ZFT samples (Z13 and Z14) with older zircon ages were located in a cooler position with respect to the ZCTA. The AFT information indicates that, in both sections, the AFT samples collected on the Present-day topography at this stage were located near the ACTA.
The frontal system uplift, in both sections, was bounded by the moderate to steep extensional faults of the Servita, Lengupa, and Tesalia complex. However, in this initial stage of the frontalsystem development, the fault offsets are not interpreted with inversion. The topographic load caused by the frontal-system initial uplift resulted in a flexural response which initiated a foreland basin (foredeep depozones, cf. Parra et al., 2009a). The development of the foreland basin accommodated the stemmed sediments proceeding primarily from the frontal-uplift erosion. The gravel-clast petrography analysis (Figure 7) presented by Parra et al. (2010) suggests that from the Middle Oligocene to ∼11 Ma, the sediments in this foreland basin came from the erosion of the Paleocene-Eocene rocks deposited in the frontal system (Barco, Los Cuervos, and Mirador Formations). Therefore, as it is represented in this model (Figure 13), the erosion window on the frontal system exposed Upper Cretaceous rocks (Guadalupe, and Chipaque Formations). This unroofing sequence and the acceptable paleoelevations for the interval between ∼23 Ma and ∼11 Ma are the most precise constraints to reconstruct the actively deforming areas for that time. The accumulative shortenings from ∼33 Ma to ∼11 Ma, are 2.4 and 4.6 km for Section 1 and 10, respectively ( Figure 13).

Restoration ca. 5 Ma
The Early Pliocene marked a depositional transition from the Lower to Upper Guayabo Formation ( Figure 14). The paleotopographic reconstruction for this model step is interpreted (for both sections) with an altitude of 2000 m in the hinterland (Wijninga, 1996;Hooghiemstra et al., 2006), which progressively decreases toward the southeast to the mean sea level in The Llanos. Different published regional interpretations inside, and nearby the study area, support the hypothesis that the deformation rates were much faster during the last 5 Ma (e.g., Martinez, 2006;Mora A. S. et al., 2013;Mora et al., 2014). In the study area  show very young AFT ages (<3 Ma) on top of a frontal basement thrust that imply acceleration of exhumation and shortening rates. To the north and along the eastern foothills cross-cutting relationships in well constrained cross-sections interpreted by Tamara et al. (2015) show that the longest traces (in profile view) of the most frontal faults cut through Late Miocene-Pliocene sedimentary units.
The analysis of the AFT samples indicates that in Section 10 the apatite ages decrease from northwest to southeast. In zone-3 and northwest part of zone-2 the apatite samples are older than 5 Ma; consequently, in our model they are located above the ACTA. From the middle of zone-2 toward The Llanos the Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 636458 samples have younger ages; therefore they are located nearby the AFTA, or even below. A similar situation is observed in Section 1 in which, apart from A35 and A36, the apatite ages from mid zone-2 toward The Llanos are younger than 5 Ma. From 11 Ma to this stage the compressional deformation continued its overall forward propagation. This caused: 1) the inversion of the master extensional faults underling the frontal system (Servita-Lengupa fault), and 2) the initiation of a series of low-angle shortcut faults, which resulted in a thin-skinned system (Boa and Guaicaramo faults). Following the early period of folding in the frontal system, in this stage the continuous shortening initiated the fault inversion of the principal extensional faults. In Section 10 ( Figure 14B), the inversion offsets for the Servita-Lengupa and Tesalia faults were 220 m and 380 m, respectively. In Section 1 ( Figure 14A), the Servita-Lengupa fault displacement was 950 m. Note that the Naranjal fault, in Section 1, has remained extensional during the entire fold-and-thrust belt formation (see frontal system in Figure 14A). A series of low angle shortcut faults broke through the stratigraphic section. They are located between the inversion faults and the foreland basin trough with an angle of approximately 30°. The seismic profiles display the low angle of the shortcut faults from the Cretaceous units to younger; however, the seismic coverage and resolution preclude a confident fault interpretation in deeper units. Therefore, these faults are tentatively interpreted here as connected to the inversion faults (cf. Rowan and Linares, 2000;Tesón et al., 2013). In Section 1, the Boa fault had 450 m of displacement, whereas the Mirador fault continued being inactive. In Section 10, a low displacement shortcut fault is interpreted at the side of the foreland basin trough; however, the most significant characteristic of zone-1 is the occurrence of the Guaicaramo fault ( Figure 14B). The Guaicaramo fault is a multi-ramp thrust fault formed by two low-angle ramps (∼30°) and a sub-horizontal segment in between, which is located in the Upper Cretaceous section. The movement of the hanging wall over the fault ramps resulted in a fault-bend fold-the Medina anticline. In this early evolutionary stage the Medina anticline (also known as Guavio anticline) displays a large wavelength and low amplitude, flanked by the Nazareth syncline to the hinterland, and the Rio Amarillo syncline to the foreland.
The gravel-clast petrography analysis (Figure 7) indicates that the erosion of the upper and middle section of the Une Formation was the main source of sediments in the foreland basin during this period . These sediments are thought to come from the uplift of the frontal system, and that fact is reflected in this evolutionary model with an erosion window exposing the lower units of the Une Formation in the proximities of the inversion faults ( Figure 14). The accumulative shortenings

Restoration ca. 3 Ma
From 5 Ma to the Late Pliocene the upper Guayabo Formation was depositing in the foreland basin ( Figure 15). In this stage the paleotopography ( Figure 6) was comparable to the Present-day elevation; with approximately 2500 m above the mean sea level in the hinterland and near the mean sea level in The Llanos (Wijninga, 1996;Hooghiemstra et al., 2006). In Section 10 ( Figure 15B), AFT samples A28 and A27 which are located in zone-1 (in the hanging wall of the inversion fault system), and A26 located in zone-2 (in the back limb of the frontal system) are younger than this period. Therefore they are located near the ACTA. Likewise, in Section 1 the AFT samples distributed from the end of the frontal-system back limb (in zone-2) to the middle of zone-1 are younger than this stage.
In Section 1 ( Figure 15A) the Servita-Lengupa inversion fault had 1450 m of displacement; whereas in Section 10 the Servita-Lengupa and Tesalia faults had 450 and 750 m, respectively. In this stage the structural activity in both sections was mostly focused in the development of the shortcut-fault system. In Section 1, the Boa fault had 4000 m of total displacement, while Mirador fault is interpreted inactive. In Section 10, the Guaicaramo fault had 6100 m of displacement and the overlying hanging wall continued its deformation increasing the Medina anticline amplitude.
The sediments deposited from 5 Ma to this stage came from the lower sections of Une Formation according to the gravel-clast petrography analysis (Figure 7). Therefore, the erosion window of this model exposes the Las Juntas Formation in the frontal system. The accumulative shortenings from ∼33 Ma to ∼3 Ma, are 9.8 and 14.5 km for Section 1 and 10, respectively ( Figure 15).

Present-Day Structural Context 0 Ma
At present, the topographic elevation in the northwest portion of both sections (in zone-3) is approximately 2800 m (Figure 16), however from the Pajarito fault toward the southeast it decreases progressively until The Llanoswhere it reaches few hundred of meters above the mean sea level (∼200 m).
During this last deformation stage, the deformation was mainly distributed along zone-1 and 2 in both sections. Folding is the dominant structural style in zone-2; whereas in zone-1 the shortcut thrusts experienced large fault slip. In Section 10 ( Figure 16B), the Guaicaramo fault accommodated 9200 m of total shortening; whereas in Section 1 the Boa and Mirador faults accommodated 8000 m and 1900 m, respectively ( Figure 16A). All the AFT samples which were located near the ACTA in zone-1 and 2 in the previous stage are at the present on the topographic surface. This indicates that during this period zone-2 and 1 underwent a rapid exhumation (approximately 4 km of exhumation in 3 Myr). This rapid exhumation is quantified and discussed below in the context of the entire deformation.

DISCUSSION
The reconstruction models highlight two marked tectonic periods: 1) post-rift subsidence during the Early and Mid Paleogene ( Figures 10A-H), and 2) positive inversion from the Oligocene to Recent (Figures 10H-R). This research emphasizes in the Oligocene to Recent contraction, which has resulted in the aforementioned tectonic inversion of this orogen.
The reconstruction models (Figures 11-16) indicate that the overall compressive deformation initiated during the Oligocene (from ∼33 to ∼23 Ma) in the hinterland (Bogota Savanna) and propagated southeastward to The Llanos. The quantitative analysis (below discussed) shows an overall constant shortening ratio of 200 m/Myr (meters per million years) from the Oligocene to Early Pliocene, which rapidly increases to approximately 2500 m/Myr from the Early Pliocene to Recent. A significant regional uplift is observed in the Bogota Savanna from the Oligovene to Early Pliocene compared to the low shortening ratio (∼200 m/Myr) during this period. This issue is discussed in further details below, as well as the non-uniqueness of this model and the structural analysis.

Non-uniqueness of the Balanced Palinspastic Restoration
Line-length balanced palinspastic restoration is a valuable tool to understand the structural evolution of a deformed geological setting; however, it offers a non-unique solution. Although the synoptic view is reasonably well constrained by the involved geological elements; any one part of the reconstructed geology may be depicted with limitations. For example, the precise structural style of individual structures cannot be established with complete certainty. In this research the kinematic evolution of the anticlines is assumed to be similar to a breakthrough (or break-thrust) anticline (cf. Willis and Willis, 1934), such is the case for the Farallones anticline in the frontal system. In this anticline an initial period of basement-involved folding from ∼23 Ma to ∼5 Ma is assumed before the extensional-fault inversion.
It has been debated whether the frontal basement uplift and antiform, the Farallones anticline, is a huge basement involved fold or an antiform whose structural relief is caused by basement duplexes (see discussion in McClay et al., 2018). The deep erosion inside the river canyons like the Guayuriba river (Mora et al., 2006) allowed to document that the basement is actually folded into a huge anticline without the presence of basement duplexes. This favors the lesser shortening solution proposed by McClay et al. (2018), which requires a distribution of less amounts of basement shortening through time.
The 2D model presented here is first order at best. It has been achieved by line-length balancing; however the internal deformation and volume loss that a thrust sheet experiences before starting to move along a discrete fault plane, in many cases, can be considerable (e.g., Williams et al., 1989). Different authors investigated the causes of volume and material loss in the crustal thicknesses of the Andes Cordillera (e.g., Eichelberger et al., 2015;Kammer et al., 2020). Eichelberger et al. (2015) ran two models using map-view reconstruction of the Andes in Bolivia to predict the crustal thicknesses of the central Andres. These models predicted volumetric excess of 8% and 25%. Eichelberger et al. (2015) argue that the volume excess in this region can be accounted for by "lower crustal flow" and/or "removal of lower crust" during the Andean deformation. Furthermore, fluid inclusions and outcrop-based structural analysis conducted in the study area  denoted that second-order folding associated to internal strain occurred prior to the first-order deformation. Kammer et al. (2020) and  document that total amounts of basement shortening by ductile deformation could represent 30% of the total shortening estimates. Additionally,  suggest that most basement folding occurred when those rocks were at depths higher than 7 km, favoring a ductile behavior of the basement. The low shortening rates from Paleogene and Early Neogene documented in this work also could be associated with this ductile deformation. However, in this research, we do not consider additional ductile or internal strain for the shortening estimates.
To constrain the kinematic evolution with thermochronometry, a constant thermal gradient of approximately 25°-30°C per kilometer in depth was assumed. This allowed for the calculation of the apatite and zircon closuretemperature approximations (ACTA and ZCTA, respectively), which are parallel to the paleotopography and located approximately 4 and 8 km beneath, respectively. However, the ACTA and ZCTA are just simplifying assumptions. The prediction of the thermochronometric ages in tectonically active settings depends on factors such as: erosion history, topographic relief, thermophysical properties of the rock, and fault kinematics and geometry. When thrust faults propagate upward along footwall ramps, the vertical component of motion and the exhumation advect heat upward. Then, whether this topographic relief is rapidly eroded, the upward advection of isotherms is enhanced. Contrarily, the sediments deposited in the foreland basin are subsequently buried to exert a downward advection of the isotherms. Therefore, assuming closure temperatures only depending on the topographic level represent a limitation of our method. Furthermore, the subsurface thermal field that is affected beneath the individual structures by both surface uplift via faulting, and erosion focused over active uplift was documented in numerous works (e.g., Ter Voorde et al., 2004;Huerta and Rodgers, 2006;Lock and Willett, 2008;McQuarrie and Ehlers, 2017). Greater accuracy of the thermal field can be obtained by means of software packages like Pecube (Braun, 2003;Whipp et al., 2009;McQuarrie and Ehlers, 2015), and Fetkin (Almendral et al., 2015;Mora et al., 2015a). Even taking all these factors into account, this evolutionary model satisfies the current known constraints of the first order structural geometry, pressure-temperature time constraint, stratigraphy, and basic shortening estimations.

Paleoelevation Data: Lipid Biomarkers Versus Paleobotanic Assessments
The scopes and challenges of using paleobotanic data to relate topographic elevations with geological ages are discussed here. As it is mentioned above, Hooghiemstra et al. (2006), and Wijninga (1996) used the nearest living-relatives-paleobotanic method to assess past elevations in the Eastern Cordillera. Different authors (e.g., Molnar and England, 1990) suggest that this is a risky method of assessing paleobotanic data, since certain species which today live at specific elevations may not have lived at the same elevations in previous times. However , Hooghiemstra et al. (2006), and Wijninga (1996) have compiled years of studies regarding the vegetation history of the Eastern Cordillera; therefore, it is difficult to debate about this knowledge without evidence that rejects the way of interpreting their data. Generally, the risk exists for the nearest-living-relatives method, but Hooghiemstra et al. (2006), and Wijninga (1996) use populations of species.
In our view, the main uncertainty of this method is not precisely that one. Instead, the risk may be to assign to the biozones from Hooghiemstra et al. (2006), and Wijninga (1996) a higher accuracy than the actual one. Perhaps this is one of the misinterpretations in , where the authors proposed that the data from Hooghiemstra et al. (2006), and Wijninga (1996) suggests a topographic growth between 6 and 3 Ma. Certainly, Figure 6, modified from the original one by Hooghiemstra et al. (2006), and Wijninga (1996), and reproduced by , does not precisely shows a topographic growth between 6 and 3 Ma, but a wider range from Mid Miocene to Pliocene.
More recently, Anderson et al. (2015) provided more accurate review about the previous hypothesis by Hooghiemstra et al. (2006), andWijninga (1996). In fact the refined ∼8 Ma age for the base of the Salto de Tequendama I and II sections by Anderson et al. (2015) is coincident with the initial assessments of Mid Miocene to Pliocene by Hooghiemstra et al. (2006), andWijninga (1996). This is also the case of the paleotemperature estimates. Hooghiemstra et al. (2006), and Wijninga (1996) do not provide a paleotemperature estimate but actually a paleoelevation range, which is broad for the Subachoque 39 and Salto del Tequendama I and II sections ( Figure 6). In such case, we can emphasize that both data sets do not disagree.

Structural Analysis and Evolution
Accumulative shortening from 33 Ma to Recent was calculated on both sections (Table 2 and Figure 17A). The shortening curves display a subparallel trend. Section 10 underwent approximately 2 km more shortening from the compression onset to Recent. Along the time, the shortening was increasing in the study area on both sections, but not constantly. From the analysis of the average accumulative-shortening curve ( Figure 17A) three marked tendencies were observed. First, from ∼33 Ma to ∼11 Ma, when the study area deformed with a relatively low and constant shortening rate of ∼160 m/Myr; experiencing an average of ∼3.5 km of accumulative shortening. Second, from ∼11 Ma to ∼5 Ma, when the rate increased to ∼650 m/Myr, resulting in an accumulative shortening of ∼7.5 km. Finally, from ∼5 Ma to Recent, when the deformation has experienced a significant acceleration that raised the rate at ∼2200 m/Myr. The total accumulated shortening at the present is 17.3 km and 19.5 km on Section 1 and Section 10, respectively. This results in an average of shortening rate of ∼550 m/Myr from ∼33 Ma to Recent.
The above described accumulative shortening was not homogeneously accommodated along the study area. Hence, to understand the propagation sequence of the deformation, as well as the strain distribution, shortening rates have been calculated for the three principal structural domains of Section 10 ( Figure 17B). These shortening rates were calculated dividing the shortening by the time frame. The following analysis is also well representative of the deformation in Section 1. The "Bogota Savanna" is one of these domains, and it covers the most northwest part of the study area to the Paja fault ( Figure 9). The "Frontal System" domain is bounded by the Paja fault to the northwest and the Tesalia fault system to the southeast. "Medina Basin" domain is flanked by Tesalia fault system to the northwest and extends to The Llanos basin.
From the beginning of the compressional deformation, approximately from 33 Ma to 23 Ma, the shortening was accommodated along the Bogota-Savanna domain. During that period, the rate of shortening was ∼260 m/Myr yielding ∼3 km of shortening (Table 3). In this period, the Frontal System and Medina Basin remained inactive. From ∼23 Ma to ∼11 Ma the structural activity was constrained in the Frontal System with a shortening rate of ∼83 m/Myr that resulted in ∼1 km of shortening; whereas the other two domains were inactive ( Figure 17B). From ∼11 Ma to ∼5 Ma the Medina Basin, which had been inactive before, rapidly started deforming with a rate of ∼1000 m/Myr. From ∼5 Ma to ∼3 Ma two domains synchronically deformed, since Medina Basin increased its shortening rate to ∼1750 m/Myr, and the Frontal System was reactivated with a rate of ∼250 m/Myr. From ∼3 Ma to Recent these two domains have remained active with rates of ∼1500 m/ Myr and ∼333 m/Myr for the Medina Basin and the Frontal System, respectively. The accumulative shortening in the Medina Basin from ∼11 Ma to Recent has been 14 km. The shortening accumulated in the Frontal System in its last reactivation (from ∼5 Ma to Recent) has been ∼2 km.
The shortening rates differentiated by structural domains clearly display an overall forward-propagation sequence of deformation in the Eastern Cordillera and Foothill. From the shortening onset (∼33 Ma) to ∼11 Ma the relatively low shortening velocities, which affected the Bogota Savanna and the Frontal System domain, are associated with the thick-skinned inversion tectonics that involved basement uplift. The structural activity from ∼11 Ma to Recent features an acceleration of the deformation. This shortening has been primarily accommodated in the thin-skinned system of the Medina Basin, and to a lesser degree in a reactivation of the Servita-Lengupa fault in the Frontal System. Parra et al. (2009a), based on a palinspastic restoration of , calculated deformation rates in the Eastern Cordillera. Like this research, they suggest the occurrence of slow and rapid (i.e., unsteadily) advance of the orogenic deformation front. Additionally, the information related to the 'orogenic migration rates' presented in Table 3 by Parra et al. (2009a) suggests that the shortening rates are in the same order of magnitude as the ones established here.

Comparison With Other 2D Palinspastic Restorations in Nearby Areas
There are two recent studies related to palinspastic restorations constrained by thermochronology near the study area: Carrillo et al. (2016) and Mora et al. (2015b). The fundamental difference  between the previous works and this investigation is that here the emphasis is on the analysis of the structural evolution, and the spatiotemporal deformation of different structural blocks; whereas the previous works focused on the refinement of the thermokinematic restorations, thermal-model validation, and testing of new technologies (Fetkin platform). Carrillo et al. (2016) restored 30 2D balanced cross-sections from Niscota (to the northeast) to Cusiana (to the southwest)located between latitude 6.30°N and 4.70°N, and longitude 73.80°W and 72.10°W. They used similar geological elements to construct the cross-sections, and because their study area is in the north, this research can be considered a continuation of that investigation. They analyzed the velocity-vector variations along geological time to understand the exhumation and thermal evolution of that area. This research presents quantitative analysis of shortening velocities and rates, which can be used to validate the Carrillo et al. (2016) hypothesis. Mora et al. (2015b), used a 2D structural section to test the capability of a new technology -Fetkin. This software gives the possibility to calibrate a kinematic restoration through the comparison between the thermal model derived from the section restoration, and the thermochronometric data collected from the field. They interpreted a larger crosssection which covers the entire Eastern Cordillera along dip direction, located south of Section 1. The thermal, exhumation, and structural evolution presented in the eastern flank of their cross-section display an overall similarity to the results of this research. They also calculated shortening rates for a number of geological periods (see Table 6 in Mora et al., 2015a) and are on the same order of the ones presented here. Note that those shortening rates were calculated in a longer cross-section (approximately the double in length), therefore approximately half of those values should be used for a comparison.
An intriguing and poorly understood aspect about the structural evolution in the study area is the high elevation and low topographic relief observed in the Bogota Savanna. These two features may indicate that a brittle deformation was not the principal mechanism that formed this highland. Similar to the two above discussed works, in our evolutionary model it was necessary to apply a sub vertical displacement to the northwest portion of the study area (i.e., Bogota Savanna). In such a way, it was possible to relate the thermal information from samples A20 and A21 to the structural evolution of this area. The structural analysis presented here suggests that the plateau formation was coincident with the initial period of slow compressive deformation from ∼33 to ∼11 Ma, in which the shortening velocity was approximately ∼160 m/Myr ( Figure 17A). Carrillo et al. (2016) dismissed the possibility of uplift by brittle deformations only such as fault-related folding (e.g.,  instead they proposed the idea of pure shear deformation associated to a mid-crust detachment. Below the brittle portion of the continental lithosphere (∼25 km) underlies the mid-lower crust (Molnar, 1988) which is characterized as being aseismic, and where the deformation may be produced by ductile processes (Allen and Allen, 2006). This ductile deformation zone has been mentioned in numerous studies associated with levels of detachments of mayor upper crustal faults (e.g., Kusznir and Park, 1987;Allen and Allen, 2006). Chamberlin (1910), proposed that the area of material in a fold that is uplifted by deformation (excess area) is equal to the product of the displacement and the depth to the detachment. Therefore, by dividing the excess area by the shortening, it is possible to obtain an estimation of the detachment depth. In this restoration, at 11 Ma the excess area below the Bogota Savanna plateau is approximately 150 km 2 , and the shortening is ∼4.6 km; therefore the estimated mid-crust detachment may be located below ∼30 km. The result of this equation may support the hypothesis of pure shear deformation proposed by Carrillo et al. (2016); however, there is a level of uncertainty in this result since the method has failed in other geological settings. The limitations of this equation have been widely discussed by numerous authors (e.g., Bucher, 1933;Kusznir and Park, 1987;Mitra and Namson, 1989;Groshong and Epard, 1994;Epard and Groshong, 1995), nevertheless, here it appears to produce a reasonable first-order result. Additionally, evidence of homogeneous flattening and strain associated with pure shear has been presented in this area previously (e.g., Kammer, 1997;Kammer and Mora, 1999;Mora and Kammer, 1999;Kammer et al., 2020). Therefore, it is possible to infer that Bogota Savanna plateau formed, during the early stage (from ∼33 Ma to ∼11 Ma), primarily by a ductile mechanism caused by the slow-deformation rate. Then, the rapid shortening velocity (from ∼11 Ma to Recent) resulted in brittle deformation primarily located in the frontal system and foreland basin, where the structural relief can be explained by the amount of shortening. Mora-Páez et al. (2016) used GPS velocities to calculate a shortening rate of ∼4000 m/Myr for the east and west portion of the Eastern Cordillera. Hence, we can estimate approximately ∼2000 m/Myr just for the eastern portion. This estimated rate is in the same order of magnitude as the ∼2200 m/Myr calculated in our sections from ∼11 Ma to Recent. However, ∼2000 m/Myr is considerably high compared to the ∼550 m/Myr calculated here from ∼33 Ma to Recent (total compressive period). Mora-Páez et al. (2016), also suggest that very young ages of deformation and low total amounts of shortening do not explain the modern GPS rates and crustal thicknesses (Poveda et al., 2015). Perhaps our shortening estimates are too low to be compared to the published estimates (Mora-Páez et al., 2016); however, if the above mentioned ductile deformation is considered, the shortening in our cross-sections could be higher. The evidences discussed here and previous researches (Mora-Páez et al., 2016) indicate that it is unlikely a Late Neogene onset of deformation in the study area.

CONCLUSION
This research proposes a comprehensive and multidisciplinary methodology to provide insights into the structural evolution of an active fold-and-thrust belt where other approaches cannot be applied (e.g., restoration guided by the growth strata). The method presented here was applied in the Eastern Cordillera of Colombia with the following outcomes.
With approximately 18 km of shortening, the study area experienced post-rift subsidence during the Paleocene -Eocene (∼65 Ma to ∼33 Ma) followed by positive inversion from the Oligocene to Recent (∼33 Ma o ∼0 Ma). During the tectonic inversion, the shortening rate increased. The compressional deformation initiated with low shortening rate (∼650 m/Myr), which increased to ∼2200 m/Myr to the Present-day. For the total amount of shortening, the deformation velocity could be considered relatively low. In particular, the low initial shortening, which has been also documented in previous studies, gives further evidence of an early ductile deformation in the study area.
The abundant evidence discussed here, regarding the timing of deformation, favors a Paleogene to Early Neogene onset of shortening and deformation along the eastern side of the Eastern Cordillera and discard a younger (Late Neogene) timings for the onset of deformation.

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