Review of geothermochronological and thermobarometric techniques for the construction of cooling and exhumation curves or paths for intrusive igneous rocks

The present study reviews radiometric and thermobarometric techniques used to construct cooling curves or paths to characterize intrusive bodies and to calculate cooling and exhumation rates. To construct these curves or paths, the temperature, time and depth variables must be estimated in intrusive bodies by applying various analytical techniques, including thermobarometry and U-Pb zircon, Ar-Ar hornblende and muscovite, fission track and (U-Th)/He zircon and apatite dating, in combination with a geological framework of reference for each intrusive body. The authors recommend to determine the crystallization age by zircon U-Pb dating, to quantify the emplacement depth using thermobarometry methods according to the composition of the intrusive body, to calculate the initial cooling ages with hornblende and muscovite Ar-Ar methods, as well as to calculate the cooling/exhumation ages in the upper crust using low-temperature thermochronology methods. These cooling curves or paths in intrusive bodies are highly relevant when studying compressive or extensional areas because they partly represent the thermal history of the era, thereby providing data on the magmatic and tectonic evolution of the region. Thus, these studies are highly important for designing geodynamic models and for their possible application in developing the tectonic model of the country.


IntroductIon
Cooling processes inside the earth's crust result from geological events, and their characterization makes it possible to learn the thermal history of a region and to obtain information about their magmatic, metamorphic and tectonic evolution (cf. Dodson, 1979;Zeitler, 1985;Green, Duddy, Laslett et al., 1989, Corrigan, 1991Gallagher et al., 1998, Reiners andBrandon, 2006). In recent decades, different radiometric dating techniques have been combined to construct cooling curves or paths in rock bodies (since isotopic systems are essentially a form of geothermometry) in order to understand the dynamics of orogenic belts and sedimentary basins, as well as their relationships with mineral deposits and petroleum systems (cf. Zeitler, 1985;Green, Duddy, Gleadow et al., 1989;Hunziker et al., 1992;Gallagher et al., 1998;McInnes et al., 1999;Armstrong, 2005;Bernet et al., 2006, Peyton andCarrapa, 2013 a, b;Bernet et al., 2019).
In regions with tectonic deformation, applying various radiometric techniques may help to constrain the exhumation time and rates, and if the exhumation results from tectonic activity, these techniques may help to determine the time of the tectonic events occurred (cf. Zeitler et al., 1982;Bernet et al., 2006;Deeken et al., 2006;Peyton and Carrapa, 2013a); however, exhumation processes are also related to erosion and climatic effects, as well as the combination of these processes to a greater extent (cf. Zeitler, 1985;Brandon et al., 1998;England and Molnar, 1990;Ring et al., 1999;Reiners and Brandon, 2006;Peyton and Carrapa, 2013a;Bernet et al., 2019).
Constructing cooling and exhumation curves for intrusive igneous rocks is an excellent opportunity to partly retrace the thermal and tectonic evolution of a region due to its extensive exposure. In addition, determining the formation age of the intrusive body can constrain deformation and post-deformation times, and offer a first-order view of the uplift or subsidence of the crust during tectonic processes Zeitler, 1985;Hurford, 1986;Anderson et al., 1988;Clarke, 1992;Anderson, 1996;Foster et al., 2001;Anderson et al., 2008;Glorie et al., 2017;Nieto Samaniego et al., 2019).
Despite the large volume of geothermochronological data on intrusive units published in the literature, few studies have focused on constructing the cooling and exhumation curves or paths of these units for understanding the evolution of areas with magmatic belts. For this purpose, the authors of this review suggest to consider that: a) reported radiometric data should be statistically validated to increase the reliability of crystallization and cooling ages of intrusive units, b) thermobarometry studies should be performed in intrusive rocks to quantify the emplacement levels of these bodies and c) carrying out systematic studies that combine the different techniques used to the construction of cooling curves or paths in intrusive rocks.
The aim of this review is to present a proposal to construct cooling curves or paths for studying the thermal history of intrusive units in Colombia, to perform a literature review on the currently used analytical techniques, and to highlight the importance of these cooling curves or paths as a complementary tool for determining the geological evolution of the respective units through geodynamic concepts. These approaches are useful to systematically calculate cooling and exhumation rates, and these inputs may be applied to complete the tectonic model of the country.

History of cooling and exhumation of intrusive igneous rocks
The exposure of large intrusive bodies in the continental crust suggests that these types of rocks play a key role in the evolution of the continents. Therefore, knowing their cooling and exhumation histories is crucial for understanding the crustal evolution that continental masses experience during tectonic events (cf. Anderson, 1996;Foster et al., 2001;Anderson et al., 2008;Mutch et al., 2016;Reiners et al., 2017, p. 197).
The history of cooling and exhumation of intrusive igneous bodies (cf. Harrison and Clarke, 1979;Hurford, 1986;England and Molnar, 1990;Foster et al., 1992;George, 1993;Caggianelli et al., 2000;Parada et al., 2005;Glorie et al., 2017;Bernet et al., 2019;Nieto-Samaniego et al., 2019) can be determined by constructing cooling curves or paths, thereby assessing the temperature, time and depth variables that characterize the body from its initial emplacement to its final exposure on the surface (Figures 1 and 2).
These cooling curves or paths are based on the fact that mineral phases dated by some isotopic method, show typical closing temperatures. Therefore, cooling ages associated with different temperature ranges can be determined in a sample, and to perform temperature vs. age plots. In addition, by thermal field modeling or by assuming a geothermal gradient in the study area, the temperature-time data can be converted into depth-time to estimate the cooling/exhumation rates that are partly associated with tectonic events (cf. Wagner et al., 1977;Zeitler et al., 1982;Zeitler, 1985;Hunziker et al., 1992;Bernet et al., 2006;Reiners et al., 2017, pp. 105-122). Figure 1. Representative diagram of a hypothetical cooling curve or path of an intrusive body from its emplacement to its final exposure on the surface, indicating the various geothermochronological and thermobarometric methods recommended for determining its cooling/exhumation history. The bars that accompany each method correspond to the uncertainty associated with age and depth. The temperature shown correspond to the hypothetical temperatures of the host rock, according to a geothermal gradient of 30 °C/km. Abbreviations: Zr: zircon; Hbl: hornblende; Ms: muscovite; Bt: biotite; ZFT: fission track in zircon; ZHE: (U-Th)/ He in zircon; AFT: fission track in apatite; and AHe: (U-Th)/He in apatite. These cooling curves or paths were initially considered empirically because they represent net histories, which involves simplifying a scenario that may have actually had more complex cooling paths (Zeitler, 1985). However, with sufficient samples and geological data, these curves or paths are highly reliable (Zeitler, 1985) using high-resolution and high-reliability ages and onsidering some type of thermal modeling (cf. Zeitler, 1985;Ehlers, 2005;Ketcham, 2005;Braun et al., 2012).
To properly construct this type of curve or path, the intrusive body should first be geologically characterized, identifying the field relationships with adjacent rocks and the type of contact observed to establish its geological context for reference ( Figure 2). This step establishes valid geological criteria that support the interpretation of the findings through various analytical methods. If the unit geological context is clear, cooling and exhumation curves or paths can be constructed by integrating various geothermochronological and thermobarometric methods and by calculating the cooling and exhumation rates of the intrusive body.

ProPosal for the constructIon of coolIng and exhuma-tIon curves or Paths of IntrusIve bodIes
The following proposal combines various geothermochronological and thermobarometric methods for determining emplacement (or crystallization) and cooling and/or exhumation ages, according to the closure temperature of each mineral system ( Figure 3). These methods also provide key data on the depth associated with each geological process. The application of these analytical methods is conditioned to the nature of the study bodies, that is, they depend on the presence of minerals, such as hornblende, biotite, muscovite, zircon, apatite, and so on. Thus, a sequence of steps is proposed below to gather the necessary data to construct the cooling curves or paths of intrusive rocks: 1. Geological characterization of the intrusive units, determining their nature, geological context and field relationships with adjacent units. 2. Geological characterization of the units adjacent to the intrusive bodies. 3. Calculation of emplacement (or crystallization) ages using the zircon U-Pb method. 4. Quantitative estimation of emplacement depths using thermobarometry techniques (e.g., aluminum in hornblende barometer), which are also applied to the host rock to estimate the geothermal paleogradient of the emplacement environment. 5. Calculation of initial cooling ages using high-temperature thermochronometers (hornblende and muscovite Ar-Ar). 6. Calculation of cooling ages in the upper crust levels using high-temperature thermochronometers (fission track and zircon and apatite (U-Th)/He dating). 7. Determination of the thermal history of intrusive bodies by combining geothermochronological and thermobarometric data, assuming a geothermal gradient or from thermal field modeling.
8. Calculation of cooling (∆temperature/∆time) and exhumation (∆temperature/geothermal gradient/∆time) rates. 9. Evaluation of the possible cooling/exhumation mechanisms associated with tectonic events, erosion or both processes.

analytIcal methods
Obtaining the necessary data to determine the cooling curve or path of intrusive rocks requires applying specific analytical methods, whose particular characteristics will influence the possible interpretations of the results. The methods used to obtain data for constructing the cooling curve or path are described below.

Zircon U-Pb geochronology: emplacement age
The emplacement age of an intrusive body can be correlated with its crystallization age, which is determined by zircon U-Pb dating (cf. Faure and Messing, 2005;Reiners et al. 2017, Figure 3. Average closure temperatures of different geothermochronometers and their potential depths, assuming a geothermal gradient of 30 °C/km. Temperature thresholds retrieved from Harrison (1982), Grove and Harrison (1996), Harrison et al., (2009) [197][198][199]. Also, the crystallization age of zircon is assumed to be the crystallization age of the igneous body because the mineral has a closure temperature of approximately 900-700 °C (Cherniak and Watson, 2001;Harrison et al., 2007) and an average of 800 °C (Figure 3), at which a melt of intermediate composition would reach 50% crystallization (Carroll and Wyllie, 1990;Harrison et al., 2007). However, the emplacement of plutonic bodies may involve multiple intrusive facies (cf. Miller et al., 2007;Reiners et al. 2017, pp. 198-199), and due to the crystallization of a zircon mineral depends on zirconium saturation and on the composition and temperature of the melt (Harrison et al., 2007) a granite body may contain zircon xenocrysts, anticrysts and autocrysts, which gives rise to a wide range of zircon U-Pb ages (cf. Coleman et al., 2004;Miller et al., 2007). Xenocrysts are zircons that are foreign to the magmatic melt, which can be derived, for example, from the source rock (in this case, termed inherited zircons) or from the host rock. Anticrysts are zircons formed in the initial magma pulses, and autocrysts are zircons that grow within the main intrusive pulse (Figure 4; Miller et al., 2007).
Based on the above information, the crystallization age of the intrusive body (or emplacement age) corresponds to the population of ages in zircons interpreted as autocrysts, with statistically reliable parameters (cf. Vermeesch, 2018). Thus, to correctly interpret the U-Pb ages in zircon crystals, the following tasks are recommended: a) characterize the mor- Xenoliths phological typology of the crystal (Pupin, 1980); b) classify internal textures using cathodoluminescence images (Corfu et al.,2003); c) cluster the ages obtained from crystals interpreted as xenocrysts, anticrysts and autocrysts; and d) calculate inherited and crystallization ages using the Isoplot (Ludwig, 2009) or IsoplotR (Vermeesch, 2018) programs.
To determine whether the calculated ages are statistically reliable, the statistical parameter mean squared weighted deviation, MSWD (Wendt and Carl, 1991), should be used because this parameter indicates the level of data dispersion (Vermeesch, 2018) and makes it possible to indirectly assess whether individual age weighting and uncertainty represent a single population (Spencer et al., 2016). Accordingly, MSWD = 1 indicates that the calculated age neither underestimates nor overestimates the associated uncertainties; MSWD < 1 indicates overestimated ages, and MSWD > 3 indicates underestimated ages (Spencer et al., 2016;Vermeesch, 2018). Thus, acceptable values of MSWD vary in the range of > 0 to 3, and these values are in turn related to the number of analyses used to calculate the age, as more analyses achieve better MSWD values (see Figure 3 of Spencer et al., 2016). For example, for a number of analyses (or spots) of N = 20, the acceptable MSWD would range from 0.35 to 1.65, for an age calculated with 2σ, which guarantees a reliable age (Spencer et al., 2016). On the other hand, statistically significant ages are those that present a value of P (X 2 ) 0.05 to 0.95, and this parameter indicates the level of data dispersion with respect to analytical uncertainties around the best-fit line (Vermeesch, 2018).

Thermobarometry in igneous rocks: emplacement depth
To approximate the emplacement depth of an intrusive body, the magma crystallization pressure and the crystallization temperature near the solidus line must be simultaneously estimated (Blundy and Holland, 1990;Holland and Blundy, 1994;Anderson, 1996;Anderson and Smith, 1995;Anderson et al., 2008;Mutch et al., 2016). Various thermobarometers have been proposed in the literature to determine these intensive variables from the composition, equilibrium or saturation of mineral phases present in the intrusive bodies (Table 1). The thermobarometers applied to intrusive rocks are described as follows:

Al-in-Hornblende barometer (Al in Hbl)
This barometer is based on the strong correlation between the total Al content in Hbl and the crystallization pressure of an intrusive (Hammarstrom and Zen, 1986). From this correlation, various equations found in the literature are formulated, generally based on experimental petrology, to estimate the crystallization pressure. The equations differ in the initial experimental parameters, such as the total composition of the host rock, the nature of the fluid phase (H 2 O or H 2 O-CO 2 ), the temperature range at which the experiment was performed and the different empirical or experimental calibrations of each equation (Hammarstrom and Zen, 1986;Schmidt, 1992;Anderson and Smith, 1995;Anderson et al., 2008;Mutch et al., 2016).
For a metaluminous granitoid of intermediate to felsic composition with a mineral association of quartz + plagioclase (An 25-35 ) + hornblende + biotite + titanium and iron oxides ± potassium feldspar ± titanite, which was most likely emplaced from 2 kbar to 14 kbar at crystallization temperatures lower than 800 °C, equation 1 (Anderson and Smith, 1995) is applied: Applying this barometer requires estimating the crystallization temperature (Anderson and Smith, 1995) because the positive correlation between temperature and [IV] Al content in Hbl (Blundy and Holland, 1990) generates overestimated values of pressure at high temperature (> 650 °C) (Anderson and Smith, 1995). To use this equation, the authors also recommend ensuring high oxygen fugacity (fO 2 ) in the system based on Fe/(Fe + Mg) ratios in Hbl ranging from 0.4 to 0.61.
For granitoids with a mineral association of amphibole + plagioclase (An 15-85 ) + biotite + quartz + alkali feldspar + ilmenite/titanite + magnetite + apatite at temperatures near the solidus line (725 ± 75 °C), which are emplaced at low pressures (from 2.5 to 0.5 kbar), equation 2 (Mutch et al., 2016) should be used: Wherein the numbers in parentheses correspond to the uncertainty 1σ expressed in the physical notation with the lowest significant number (Mutch et al., 2016). To correctly apply this equation, the crystallization temperature should also be estimated, thereby ensuring the necessary conditions for its use, which assumes a crystallization temperature of approximately 725 ± 75 °C, as previously mentioned.

Amph+Pl Thermometer
The Amph + Pl thermometer is used to determine the crystallization temperature of a granitoid of variable composition (even with rocks unsaturated in quartz) with a mineral association between amphibole and plagioclase over a temperature range of 400-1000 °C and a pressure range of 1-15 kbar (Holland and Blundy, 1994). This thermometer is based on the equations of equilibrium between edenite + quartz = tremolite + albite (for associations with quartz, termed thermometer A) and edenite + albite = richterite + anorthite (for associations without quartz, termed thermometer B). This technique considers a non-ideal mixture in both amphibole and plagioclase compositions and calibrates the equations using a database of natural and synthetic amphiboles to reach uncertainties of ± 35-40 °C or higher for iron-rich amphiboles (Holland and Blundy, 1994). The equations (3) and (4) are proposed by the authors.
Iterative calculations using this thermometer and the Al-in-Hbl barometer can provide reliable pressure and temperature ranges in metaluminous granitoids (cf. Anderson et al., 2008).
Thermobarometers based on the composition of calcic amphiboles A set of barometers and thermometers based on the chemical composition of amphiboles, which is applicable under conditions of temperature raging from 1130 to 800 °C and under a wide pressure range from 22 to 1.3 kbar, have been reported for calc-alkaline and alkaline rocks with calcic amphibole (Ridolfi et al., 2010;Ridolfi and Renzuli, 2012). These thermobarometers are empirically formulated from data in the literature regarding natural and synthetic amphibole compositions, with multivariate least squares analysis of experimental compositions of this mineral (Ridolfi et al., 2010;Ridolfi and Renzuli, 2012).
Thermobarometers based on the clinopyroxene-liquid equilibrium Thermobarometers based on the clinopyroxene-liquid equilibrium can be used in intrusive rocks of mafic composition under conditions of pressure ranging from 0 to 30 kbar and temperatures ranging from 1100 to 1475 °C (Putirka et al., 1996). These thermobarometers are recommended for rocks that are not granitoids, because mafic rocks are exposed over wide areas in Western Colombia.

Other barometers
In peraluminous granitoids without the required mineral associations to use the Al-in-Hbl barometer and without a suitable composition of calcic amphiboles, pressure can be estimated using empirical garnet + biotite + plagioclase + quartz (Wu et al., 2004a) or garnet + muscovite + plagioclase + quartz (Wu et al., 2004b) barometers developed for metapelitic rocks, which are recommended for granitoids with this mineralogy (Anderson et al., 2008). Similarly, in peraluminous plutons with garnet and muscovite, where muscovite contains celadonite, pressure can be estimated by barometry developed in high-pressure metapelites (cf. Wei and Powell, 2004;2006;Anderson et al., 2008).
Moreover, the identification of stable mineral phases under specific pressure conditions can constrain the calculation, such Thermometer A: ( 3) where the term Y ab = 0, if > 0.5; otherwise, Y ab = 12.0(1 − X ab ) 2 − 3.0 kj; the empty box in X A corresponds to the vacancy in the crystallographic site A of amphiboles.
Thermometer B: where the term Y ab−an = 0, if X ab > 0.5; otherwise, Y ab−an = 12.0 (2X ab − 1) + 3.0 kj, expressing the temperature as Kelvin, the pressure as kbar, the composition X i Φ as the molar fraction of element i in mineral phase Φ , and gas constant R = 0.0083144 kjK −1 mol −1 (Holland and Blundy, 1994 as the presence of magmatic epidote, which suggests a pressure > 5 kbar (Zen and Hammarstrom, 1984;Schmidt, 1993). The emplacement depth of an intrusive body can also be estimated when knowing the pressure and temperature conditions at the contact aureole (e.g., garnet + biotite + cordierite equilibrium [cf. Treloar, 1981] and with the aforementioned barometry in peraluminous granitoids applied to the contact rocks [cf. Wu et al., 2004 a, b]). Moreover, the depth of an intrusive body can be qualitatively estimated in terms of epizonal (0-10 km), mesozonal (6-16 km) and catazonal (> 10 km) emplacement by characterizing the type of contact and by describing microand mesostructures of the intrusive body and host rock (Buddington, 1959;Paterson et al., 1996), although the resulting ranges should be taken with caution.

Other thermometers
Other thermometers, which can be applied to intrusive bodies of different compositions, are based on the saturation of accessory phases or on the content of the trace elements of these phases (cf. Anderson et al., 2008), such as zircon (Watson and Harrison, 1983) and apatite (Tollari et al., 2006) saturation thermometers and titanium (Ti)-in-zircon (Watson et al., 2006;Ferry and Watson, 2007) and zirconium (Zr)-in-titanite (Hayden et al., 2008) thermometers. Crystallization temperatures obtained by zircon saturation and Ti-in-zircon thermometers are not always comparable, and the crystallization temperatures tend to be more closely associated with the temperature derived from the Ti-in-zircon method (Harrison et al., 2007) because zircon saturation temperature can be affected by the presence of pre-magmatic zircons (xenocrysts) in the granitic systems that saturate the magma with zircon during the melting process and not during the crystallization process.
The Ti-in-zircon thermometer is based on both empirical and experimental calibration, according to the abundance ratio of titanium in zircon and to its variation as a function of the crystallization temperature, and it is derived from equation 5 (Watson et al., 2006): However, the substitution of Ti in zircon is related to the activities of TiO 2 and SiO 2 . Therefore, a new experimental calibration is required in a temperature range from 1450 to 590 °C, which emphasizes the effect of these activities (Ferry and Watson, 2007): (6) In turn, the Zr-in-titanite thermometer was empirically and experimentally calibrated by Hayden et al. (2008), based on experiments at temperatures ranging from 800 to 1000 °C and at pressures ranging from 10 to 24 kbar, and the variation in the Zr abundance in titanite as a function of temperature and pressure is described by using equation 7:  Watson and Harrison (1983) Apatite-saturation thermometer Near-liquidus temperature Peraluminous granitoids Tollari et al. (2006) Ti-in-zircon thermometer Temperature range = 1450-590 °C Intermediate to felsic granitoids Ferry and Watson (2007) Zr-in-titanite thermometer Pressure range = 10-24 kbar, temperature range = 800-1000 °C Intermediate to felsic granitoids Hayden et al. (2008) This thermometer may underestimate the temperature if the TiO 2 content is < 1, as expected in common intrusive bodies without rutile (Hayden and Watson, 2007;Anderson et al., 2008). Also, the apatite saturation temperature can be used for peraluminous granitoids (Tollari et al., 2006).

High-temperature thermochronology: initial
cooling age The interpretation of a cooling age in intrusive rocks is directly related to the concept of closure temperature. Closure temperature is defined as the temperature at which a radiometric system starts to accumulate the daughter isotopes produced by the decay of the parent isotope, which reduces the diffusion loss to almost zero (Dodson, 1973;1979). Based on this concept, various thermochronometers are sensitive to temperature variations and provide information on the thermal history of the rock, rather than mineral crystallization ages (cf. Dodson, 1973;Zeitler, 1985;Gallagher et al., 1998;Kelley, 2002;Reiners and Brandon, 2006;Peyton and Carrapa, 2013a).
Thus, the cooling ages of intrusive rocks correspond to the ages obtained using thermochronological systems with a closure temperature lower than the crystallization temperature of the intrusive body (cf. Dodson, 1979;Zeitler, 1985;Kelley, 2002;Reiners and Brandon, 2006). The initial cooling ages correspond to the post-emplacement ages that are close to the crystallization age of the intrusive body. To calculate these ages, high-temperature thermochronology is the most suitable approach (cf. Harrison, 1982;Zeitler, 1985;Reiners and Brandon, 2006;Schaen et al., 2020). To clarify, the closure temperatures of the systems will depend on the cooling rate of the rock and on kinetics, geometry, composition and size crystal (cf. Dodson, 1979;Harrison, 1982;Donelick et al., 2005;Reiners and Brandon, 2006;Reiners et al., 2017, pp. 97-101;Oriolo et al., 2018;Bernet et al., 2019).
Among the high-temperature thermochronometers (Figure 3), hornblende and muscovite 40 Ar/ 39 Ar methods stand out because the radiogenic isotope argon is a noble gas, which can be easily affected by diffusion in the minerals and which is only fully retained when the intrusive rock is already at a cooling stage (Dalrymple and Lanphere, 1969). Furthermore, these thermochronometers have considerably higher closure temperatures than other thermochronometers suitable for calculating initial cooling ages (Harrison, 1982;Zeitler, 1985;Harrison et al., 2009;Reiners et al., 2017, pp. 231-253;Schaen et al., 2020) (Figure 3). In the biotite 40 Ar/ 39 Ar method, the thermometer is considered to be a low-temperature thermometer (Figure 3) because the closure temperatures range from 330 to 280 °C (cf. Harrison et al., 1985;Grove and Harrison, 1996), thus recording younger cooling ages, which would not correspond to the initial cooling ages of an intrusive rock. The potassium feldspar 40 Ar/ 39 Ar method shows a wide range of closure temperatures (~350-150 °C) ( Figure 3) and frequently contains complex microstructures and degrees of alteration. Thus, the potassium feldspar 40 Ar/ 39 Ar thermochronometer is not recommended for plutonic and metamorphic rocks and is accordingly more often applied to the geochronology of volcanic rocks (cf. Reiners et al., 2017, p. 240;Schaen et al., 2020). 40 Ar/ 39 Ar dating can lead to interpretation problems and ages without geological significance. These caveats are related to a) 39 Ar recoil, which disturbs the 40 Ar/ 39 Ar ratio and the accuracy of age and affects fine-grained minerals and minerals with a high surface/volume ratio (e.g., micas) (cf. Kelley, 2002;Reiners et al., 2017, p. 240;Schaen et al., 2020); b) low K concentrations (~< 2% wt K 2 O), which makes it difficult to detect 39 Ar, and high Ca/K ratios, which can generate interference problems due to reactions in the reactor, two conditions that affect amphiboles (cf. McDougall and Harrison, 1999;Kelley, 2002;Reiners et al., 2017, p. 241); c) argon loss by diffusion, due to the density of defects in the crystalline structure of the mineral (e.g., fractures) and to edge defects (cf. Dahl, 1996;Kramar et al., 2001;Kelley, 2002;Schaen et al., 2020), and d) argon excess in the minerals resulting from fluid, melt or solid inclusions, with the fluid and melt phases being high reservoirs of 40 Ar because the element is highly incompatible (cf. Esser et al., 1997;Boven et al., 2001;Kelley, 2002;Schaen et al., 2020).
Considering the above information, step heating should be used in 40 Ar/ 39 Ar dating because this approach may show the behavior of a partly open system, the presence of non-radiogenic 40 Ar in the sample and the 39 Ar recoil associated with the surface of the crystal (cf. Turner, 1970;Kelley, 2002;Reiners et al., 2017, pp. 243-246;Schaen et al., 2020). This technique generates an age spectrum as a function of the degassed fraction of the sample from step heating, in which a spectrum with a significant number of continuous steps and with indistinguishable ages will provide a plateau age, and this age is considered to be a concordant age (cf. Lanphere and Dalrymple, 1971;McDougall and Harrison, 1999;Fleck et al., 1977;Faure and Mensing, 2005;Reiners et al., 2017, pp. 243-246;Schaen et al., 2020). The plateau age should comprise at least 50% of the 39 Ar released, at a 95% confidence level for each age encompassing the plateau age (Fleck et al., 1977), as well as include at least five or more consecutive steps (Schaen et al., 2020).
The spectrum of ages shows the excess of 40 Ar in the samples, for example, the excess of 40 Ar accumulated through fluid inclusions when released at low temperatures, as well as the excess of 40 Ar accumulated through melt or solid inclusions when released at high temperatures, which generates a U-or horseshoe-shaped spectrum of ages with older ages at the ends (cf. Lanphere and Dalrymple, 1971;Wartho et al., 1996;Kelley, 2002). Spectra with discordant ages that do not form plateau ages are also observed and are attributed to 40 Ar loss by diffusion, 39 Ar recoil or 39 Ar excess (cf. Turner, 1970;Boven et al., 2001;Reiners et al., 2017, pp. 243-246;Schaen et al., 2020).

Initial cooling rate
The initial cooling rates can be calculated using the hornblende and muscovite 40 Ar/ 39 Ar method thanks to the different closure temperatures of these thermochronometers (Figure 3). In this case, the oldest cooling age would be recorded when applying the hornblende 40 Ar/ 39 Ar method, which has a higher closure temperature (between 550 and 480 °C) (cf. Harrison, 1982) than the muscovite 40 Ar/ 39 Ar method, which has a relatively lower closure temperature (between 440 and 390 °C) (cf. Harrison et al., 2009), thereby yielding younger ages. By calculating these ages, including the crystallization age, the initial cooling rates of the intrusive rock can be estimated.
This initial cooling rate can be interpreted as a) cooling associated with heat transfer by conduction or advection (cf. Peacock 1989;Zeitler, 1985;Ehlers, 2005), or b) cooling attributed to tectonic exhumation processes that decrease the ambient temperature of the intrusive body (cf. Zeitler, 1985;England and Molnar, 1990;Ring et al., 1999). To infer the process that controlled the initial cooling, the tectonic context in which the intrusive body was emplaced and the physical parameters that describe the heat transfer processes should be considered (Ehlers, 2005). For example, one-dimensional heat diffusion in a given intrusive body could be estimated using equation 8 (Carslaw and Jaeger, 1959;Ehlers, 2005): where ∂T is the temperature difference between the intrusive body and the host rock, ∂t is the time that the system spends decreasing its temperature, z is the parameter associated with the dimensions of the intrusive body, and K is the thermal diffusivity constant of an intrusive body (approximately 32 km 2 / Ma) (Ehlers, 2005). The dimensions of the intrusive body can be estimated with its exposed surface area, field relationships and depth of emplacement, which is combined with geophysical techniques if available (cf. Clarke, 1992;Pitcher, 1997). The variation in temperature as a function of time can be assessed using the previously described thermometers and thermochronometers, including thermobarometry in the host rock, to determine the possible thermal paleogradient. By applying equation 8 to a one-dimensional model and assuming a flat body shape, the initial temperature of a surface intrusive body with a diameter of 5 km, which affects a host rock at a temperature of 50 °C, decreases from 700 to 500 °C over 100 000 years, which is exclusively due to heat conduction processes (Ehlers, 2005). In this model, the intrusive body size (diameter) increases the amount of heat that affects the host rock as its temperature decreases, as well as the duration and distance of the heat source. Under these conditions, an equilibrium should be reached between the temperature of the intrusive body and the host rock, which is controlled by the geothermal gradient of the environment (Ehlers, 2005).
Considering the above information, simple heat transfer models should be applied to identify the prevailing cooling process in each intrusive body, that is, whether cooling results from heat dissipation, from an associated tectonic event or from both (cf. Zeitler, 1985;Ehlers, 2005). For this, it must determine the cooling rates, the closure temperatures of the geothermochronometers, the initial temperatures of both the intrusive body and the host rock (thermobarometry, Table  1) and the relative dimensions of the rock bodies (cf. Clarke, 1992;Pitcher, 1997;Ehlers, 2005).

Low-temperature thermochronology: cooling/ exhumation age in the upper crust
Low-temperature thermochronology is suitable for studying cooling processes at upper levels in the crust (cf. Gallagher et al., 1998;Farley, 2002;Peyton and Carrapa, 2013a;Reiners et al., 2017, pp. 105-127;Bernet et al., 2019) because this approach includes thermochronometers with closure temperatures lower than 300 °C, equivalent to shallow depths not exceeding 10 km, for a geothermal gradient of 30 °C/km (Figure 3, Table 2). In general, this cooling is generated as rocks move toward the surface through exhumation events resulting from erosion processes or tectonic activity (cf. England and Molnar, 1990;Brandon et al., 1998;Ring et al, 1999;Reiners and Brandon, 2006;Peyton and Carrapa, 2013a). These ages can be associated with the exhumation ages of the rock body. By combining different thermochronometers (with different closure temperatures), exhumation rates can be estimated and thermal histories can be modeled during the passage of rocky bodies through the upper crust (Reiners and Brandon, 2006;Peyton and Carrapa, 2013a;Reiners et al., 2017, pp. 105-127;Bernet et al., 2019;Nieto-Samaniego et al., 2019).
The most commonly used low-temperature thermochronology techniques (cf. Gallagher et al., 1998;Farley, 2002;Donelick et al., 2005;Reiners and Brandon, 2006;Chew and Spikings, 2015;Bernet et al., 2019) are zircon and apatite (U-Th)/ He (ZHe and AHe) and fission track (ZFT and AFT) dating. To correctly apply these thermochronological methods, the concepts of partial retention zone (PRZ), applicable to the (U-Th)/He technique, and partial annealing zone (PAZ), applicable to the fission track technique (Gleadow and Fitzgerald, 1987;Wolf et al., 1998;Brandon et al., 1998;Armstrong, 2005;Reiners and Brandon, 2006;Peyton and Carrapa et al., 2013a;Bernet et al., 2019) should be clarified, in addition to the concept of closure temperature. These additional concepts must be considered because the mineral may experience different cooling processes, including a) constant, monotonic and relatively fast cooling with a temperature that is always decreasing over time, in which case the concept of closure temperature would be applied, thereby assessing the cooling age (Table 2, Figure 5A), and b) variable cooling over time, possibly remaining in the temperature ranges of the PAZ and PRZ, in which case the systems are partly open, thus generating partial reset ages (Table 2, Figure 5B) (cf. Gallagher et al., 1998;Wolf et al.,1998;Brandon et al., 1998;Farley, 2002;Donelick et al., 2005;Reiners and Brandon, 2006;Bernet et al., 2019). Hence, thermochronological systems can follow different cooling paths, and each technique must be adequately applied toward understanding their thermal history and modelling cooling and exhumation rates close to reality (cf. Gallagher et al., 1998;Wolf et al., 1998;Farley, 2002;Reiners and Brandon, 2006;Bernet et al., 2019).
(U-Th)/He dating is based on the accumulation of He particles in zircon and apatite crystals derived from alpha decay in 238 U → 206 Pb, 235 U → 207 Pb, and 232 Th → 208 Pb decay systems (Wolf et al., 1996). Because He undergoes diffusion processes at high temperatures, He is only retained in the crystal below the closure temperature of each mineral (Figure 3, Table 2), thereby recording the cooling age (Wolf et al., 1996;Farley, 2000;Farley, 2002), or He might be partly accumulated in the partial retention zones (PRZ) established for each thermochronometer (Table 2), whereby the ages would be considered to be partial reset ages (Wolf et al., 1998).
Fission track dating is based on the spontaneous fission of 238 U that occurred in the crystalline structure of the minerals that contain U (Fleischer et al., 1975). This fission generates tracks (or damage) in the crystalline structure that are formed constantly over time, as long as they are generated below the closure temperature of each system, because fission tracks can quickly annealed due to solid-state diffusion in high-temperature environments (cf. Gallagher et al., 1998;Donelick et al., 2005;Reiners and Brandon, 2006;Bernet et al., 2019). The density of tracks in a crystal sector and the concentration of 238 U make it possible to calculate the age at which the mineral cooled below the closure temperature, according to a thermochronometer (Table 2), as long as the mineral has experienced relatively fast, monotonic and constant cooling at rates of approximately 10 to 100 °C/Ma (cf. Bernet et al., 2019) ( Figure  5A). Conversely, if the minerals experienced slow cooling (2 °C/Ma), staying in the range of PAZ temperatures (Table 2, Figure 5b), or overheating due to a thermal event (in intrusive bodies), the tracks can partially annealed, thus recording partial reset ages (cf. Bernet et al., 2019). The AFT technique, in addition to calculating ages, requires studying/understanding of partial annealing degree experienced by the fission track (that is, how much they are shortened or anneling) necessary to know if the sample remained at PAZ temperatures, and for modeling the thermal history of the mineral (Gleadow et al.,1983;Green et al., 1986;Green, Duddy, Laslett et al., 1989;Gallagher et al., 1998;Ketcham et al., 1999;Donelick et al., 2005;Peyton and Carrapa, 2013a). The partial annealing process is mainly controlled by the time-temperature factor that affected the sample and by the direction of tracks with respect to the crystallographic axis, which is correlated with apatite composition (Gleadow et al., 1983;Green et al., 1986;Ketcham et al., 1999;Donelick et al., 2005). The temperature-time factor is determined by the distribution of horizontally confined track lengths, and the crystallographic direction of tracks (and apatite composition) is obtained using kinetic parameters (Gleadow et al., 1983;Green et al., 1986;Green, Duddy, Laslett et al., 1989;Gallagher et al., 1998;Ketcham et al., 1999;Donelick et al., 2005;Peyton and Carrapa, 2013a;Bernet et al., 2019).
The track length distribution method is based on the length-temperature relationship that exists in fission tracks, which could be shortened if the mineral resides within the PAZ temperatures for a significant time, consequently showing a bimodal track length distribution (cf. Green et al., 1986;Donelick et al., 2005;Peyton and Carrapa, 2013a). In this bimodal distribution, the shortest tracks would be the oldest, and the longest tracks would be the most recent, which are formed below the PAZ temperatures (cf. Green et al., 1986;Donelick et al., 2005;Peyton and Carrapa, 2013a).
One of the kinetic parameters is the chemical composition of apatite, which is based on a strong track annealing dependence on the chlorine/fluorine ratio, whereby chlorine-rich apatites are more resistant to track annealing, whereas apatites with fluoride experience faster partial annealing (Green et al., 1985;Barbarand et al., 2003;Donelick et al., 2005). Because determination of the apatite composition requires using destructive or complex and costly analytic techniques, an alternative kinematic parameter is Dpar (Donelick, 1993;Barbarand et al., 2003;Donelick et al., 2005).
Dpar is based on the solubility of the crystal to acid digestion to reveal tracks (etching) and is defined as the diameter of the geometric figure formed by the intersection between the fission track and the polished surface of the crystal; low Dpar values are typical of fluorine apatites, characterized by low annealing temperatures, whereas high Dpar values are common-a b Figure 5. Possible cooling/exhumation curves of paths in the upper crust. a) Relatively fast and monotonic cooling curve; age derived from the closure temperature step (cooling age calculation). b) Slow cooling curve with prolonged passage in the PAZ and PRZ (partial reset age calculation). Temperature thresholds retrieved from Bernet et al. (2019) and references therein (see Table 2). ly found in chlorine-rich apatites, characterized by high annealing temperatures (Donelick, 1993;Barbarand et al., 2003;Donelick et al., 2005). In summary, the distribution of horizontally confined fission track lengths and the kinetic parameter Dpar should be determined to modeling the thermal history of apatite (cf. Green et al., 1985;Donelick, 1993;Gallagher et al., 1998;Barbarand, 2003;Donelick et al., 2005).
In the ZFT method, partial track annealing is controlled by the temperature, time and accumulation of alpha radiation damage (cf. Nasdala et al., 2001;Rahn et al., 2004;Reiners and Brandon, 2006). Alpha radiation damage occurs during 238 U, 232 Th and 235 U decay and causes changes in the crystallography of the mineral (or metamictisation), as well as color accumulation in zircon due to point defects in the crystals (Nasdala et al., 2001;Kamp, 2002, Garver et al., 2002;Rahn et al., 2004).
Various studies suggest that zircons with alpha radiation damage are more likely to show partial annealing at temperatures below those faced by zircons with little radiation damage and lower closure temperatures (cf. Brandon et al., 1998;Rahn et al., 2004;Reiners and Brandon, 2006). Thus, experiments with zircons without alpha radiation damage indicate a closure temperature of approximately 342 °C, at a cooling rate of 10 °C/Ma (Rahn et al., 2004;Reiners and Brandon, 2006), which is considerably higher than that estimated in zircons with natural radiation damage of approximately 230 °C under the same cooling conditions (Brandon et al., 1998) (Table 2).
The closure temperature of zircons with radiation damage is the most recommended value for most geological processes (Brandon et al., 1998;Reiners and Brandon, 2006). In turn, the closure temperature of zircons with zero radiation damage can be applied to environments where zircons show rapid cooling and derive from high-temperature environments (> 350-400 °C), causing little alpha radiation damage to these zircons during closure (Reiners and Brandon, 2006).

Cooling and partial reset age
Various cooling ages can be calculated when an intrusive body experiences a monotonic and relatively fast cooling ( Figure  5a). Based on the age vs. closure temperature correlation, the oldest ages would be recorded by ZFT and ZHe (closure temperatures that would vary between approximately 230 and 190 °C, respectively; Table 2), which would provide data on the time when the intrusive passed to a depth ranging from 7 to 5.6 km, assuming a geothermal gradient of 30 °C/km and a surface temperature of 20 °C (cf. Gallagher et al., 1998;Farley, 2002;Reiners and Brandon, 2006;Peyton and Carrapa, 2013a;Bernet et al., 2019). Similarly, the youngest ages would be assessed with AFT and AHe (Table 2), which would indicate the time when the rock body passed at a depth of 3 km with AFT and 1.6 km with AHe (depths close to the surface) (cf. Gallagher et al., 1998;Farley, 2002;Reiners and Brandon, 2006;Peyton and Carrapa, 2013a;Bernet et al., 2019). Therefore, when combining the four thermochronological techniques (or at least two), the cooling ages associated with exhumation processes can be calculated at upper crust levels, as well as the cooling/exhumation rates (cf. Zeitler et al., 1982;Zeitler, 1985;Brandon et al., 1998;Ring et al., 1999;Reiners and Brandon, 2006;Deeken et al., 2006).
Conversely, if the rock body experienced slow and variable cooling ( Figure 5) or reheating events, such as through magmatism, hydrothermalism or burial metamorphism, it could reside in the PAZ and PRZ, and the ages calculated using the ZFT, ZHe, AZT and AHe methods would correspond to partial reset ages, whose interpretation requires direct or indirect thermal modeling based on numerical methods (Ehlers, 2005;Peyton and Carrapa, 2013a;Bernet et al., 2019). The literature contains programs and routines for modelling the thermal history of one or more samples based on thermochronological data, such as the HeFty (Ketcham, 2005), QTQt (Gallagher, 2012) and Pecube (Braun et al., 2012) programs. The last program solves the heat transfer equation in three dimensions, which reduces exhumation scenarios, tectonic configurations and time-temperature paths (thermal history) of study samples (cf. Braun et al., 2012).
Considering the above information, two or more low-temperature thermochronometers should be applied when studying the cooling/exhumation processes that intrusive bodies have experienced in their passage through the upper crust, and to thermal modeling when the calculated ages are interpreted as a product of partial reset.

conclusIons and fInal consIderatIons
The path of an intrusive body from emplacement to final exhumation at the surface should be studied to determine its cooling and exhumation history. For this purpose, geothermochronology and thermobarometry techniques can be combined with field observations and relationships.
The present review highlights the application of geothermochronological and thermobarometric methods to determine temperature-time-depth variables in intrusive bodies.
Granitoid emplacement ages are determined using the zircon U-Pb method, and their depths can be estimated with thermobarometry techniques (Table 1). Initial cooling ages can be calculated with high-temperature thermochronology techniques, including hornblende and muscovite 40 Ar/ 39 Ar methods. Initial cooling rates can also be calculated by combining the crystallization age with the initial cooling age. These cooling rates can be associated with heat diffusion processes (by conduction or advection) or with tectonic events.
Low-temperature thermochronology techniques can be applied to study the passage of intrusive bodies in the upper crust level (<10 km): zircon and apatite fission track (ZFT and AFT) and (U-Th)/He (ZHe and AHe) dating. The ZFT and ZHe thermometers provide data on cooling/exhumation ages at depths ranging approximately from 8 to 6 km (assuming a geothermal gradient of 30 °C/km), and AFT and AHe thermometers record cooling/exhumation ages at near-surface depths (approximately from 4 to 2 km). If the mineral has remained in the PAZ and/or PRZ zones, the calculated ages are considered to be partial reset ages and require thermal modeling for correct interpretation. Two or more low-temperature thermochronology techniques should be combined to calculate exhumation rates at upper crust levels and their possible correlations with tectonic events.
The different techniques proposed in this review may present sources of uncertainty that are not normally included in the error propagation calculation, such as sampling bias, number of statistically representative data, human errors in laboratory work and interlaboratory standard calibration errors, which increase the uncertainty around the final curve or path.
Finally, establishing the cooling/exhumation history of an intrusive body is crucial for understanding the thermal and tectonic evolution of a region. Furthermore, determining the cooling and exhumation rates helps to limit the main tectonic pulses associated with the construction of an orogenic belt or with its collapse during crustal evolution.

acknowledgements
This study was conducted by the Tectonics Group of the Technical Dirección técnica de Geociencias Básicas of the Servicio Geológico Colombiano (SGC) within the framework of the project entitled Modelo Tectónico de Colombia [Tectonic Model of Colombia]. We would like to sincerely thank the anonymous reviewers for their valuable comments and suggestions, which helped to improve the manuscript. We also thank geologists Carlos Augusto Quiroz Prada and Ana María Patiño for their important suggestions and technical comments.