Keywords: Remote sensing; Thermal infrared; Temperature; Geothermal; Prospects; Paipa. Palabras clave: Sensores Remotos; Infrarrojo Termal; Temperatura; Geotermia; Prospectos; Paipa. How to cite item Calderon-Chaparro, R. A., & Vargas-Cuervo, G. (2019). Determination of Hydrothermal Prospects in the Geothermal Region of Paipa (Boyacá, Colombia), Using Remote Sensing and Field Data. Earth Sciences Research Journal, 23(4), 265-282. DOI: https://doi.org/10.15446/ esrj.v23n4.77810 Geothermal resources have been found with the help of conventional geological, geochemical, and geophysical techniques. However, in some cases, difficult access to research areas hinders its application, which in turn increases operating costs and restricts areas of exploration. Remote sensors with thermal infrared are an important tool in the exploration of geothermal resources due to the low cost and high efficiency in the study of large geographic areas. In this research, thermal images of remote satellite sensors are used and integrated with existing geological-geophysical information for the determination of exploratory prospects of hydrothermal vents in the geothermal region of Paipa, Boyacá. The satellite thermal images used are Landsat 7 ETM+ and Landsat 8 OLI/TIRS. Radiometric calibration, atmospheric correction, surface thermal emissivity calculation, a single-channel, and Split window algorithm were applied to retrieve the land surface temperature. Soil temperature, surface temperature, air temperature, relative humidity, atmospheric pressure, and thermal imagery of the geothermal manifestations were also used. The field data contributed to the correction of the thermal bands, and the soil temperature data helped to create a subsurface temperature map at one-meter depth. Once the primary and secondary data was available, in a GIS, an unweighted spatial model was implemented, which used four input indicators (satellite temperature index, soil temperature index, structural lineaments index, and iso-resistivity index), to determine the areas most likely to find geothermal fluids. Six areas of exploratory interest for the extraction of hydrothermal fluids were determined, of which two are already known, and the remaining four are proposed as prospects. The results allowed us to conclude that thermal remote sensors are a useful tool for the mapping of geothermal anomalies in the Paipa region, and by combining these anomalies with geological-geophysical information, it is possible to determine specific areas of exploration. ABSTRACT Determination of Hydrothermal Prospects in the Geothermal Region of Paipa (Boyacá, Colombia), Using Remote Sensing and Field Data Determinación de Prospectos Hidrotermales en la Región Geotérmica de Paipa (Boyacá, Colombia), Utilizando Sensores Remotos y Datos de Campo ISSN 1794-6190 e-ISSN 2339-3459 https://doi.org/10.15446/esrj.v23n4.77810 Los recursos geotérmicos han sido encontrados con ayuda de técnicas de campo, geológicas, geoquímicas y geofísicas, las cuales, en algunas ocasiones, el difícil acceso al área de investigación, dificulta su aplicación, aumentando los costos de operación y restringiendo espacialmente las áreas de exploración. Los sensores remotos con infrarrojo termal, son una herramienta importante en la exploración de recursos geotérmicos, debido al bajo costo y alta eficiencia en el estudio de grandes áreas geográficas. En esta investigación se utilizaron e integraron imágenes termales de sensores remotos satelitales, con información geológica-geofísica existente, para la determinación de prospectos exploratorios de fuentes hidrotermales en la región geotérmica de Paipa, Boyacá. Las imágenes térmicas satelitales utilizadas son Landsat 7 ETM+ y Landsat 8 OLI/TIRS. También, se utilizaron datos de campo, de temperatura de suelo, temperatura de superficie, temperatura de aire, humedad relativa, presión atmosférica e imágenes térmicas de las manifestaciones geotérmicas. Las bandas termales Landsat fueron calibradas radiométricamente, posteriormente corregidas por efectos atmosféricos y emisividad de la superficie, aplicando algoritmos de un solo canal y Split Window, para Landsat 7 ETM+ y Landsat 8 TIRS, respectivamente. Los datos de campo contribuyeron a la corrección de las bandas termales, y los datos de temperatura de suelo, ayudaron a la creación de un mapa de temperatura subsuperficial a un metro de profundidad. Una vez se contaba con los insumos primarios y secundarios, en un SIG, se implementó un modelo espacial no ponderado, el cual utilizó cuatro variables de entrada (índice temperatura satelital, índice temperatura de suelo, índice lineamientos estructurales e índice iso-resistividad), para determinar las zonas de mayor probabilidad de encontrar fluidos geotérmicos. Se determinaron seis áreas de interés exploratorio para la extracción de fluidos hidrotermales, de las cuales, dos son ya conocidas y las cuatro restantes se proponen como prospectos. Los resultados permitieron concluir que los sensores remotos termales son una herramienta útil para el mapeo de anomalías geotérmicas en la región de Paipa, y combinando estas anomalías con información geológica-geofísica, es posible determinar áreas puntuales de exploración. RESUMEN Record Manuscript received: 11/02/2019 Accepted for publication: 02/08/2019 EARTH SCIENCES RESEARCH JOURNAL Earth Sci. Res. J. Vol. 23, No. 4 (December, 2019): 265-282 Rafael Andrés Calderón-Chaparro, German Vargas-Cuervo Universidad Nacional de Colombia * Corresponding author: raacalderonch@unal.edu.co G EO LO G Y 266 Rafael Andrés Calderón-Chaparro, German Vargas-Cuervo Introduction Geothermal energy refers to the heat produced inside the Earth. This energy is transported to the surface through water or steam (DiPippo, 2015), and is considered one of the cleanest and most renewable energy sources on the planet (Dickson & Fanelli, 2003). Geothermal fluids are used to produce electricity and have several applications, including space heating, geotourism, and health, among others (Lund & Boyd, 2015). The geothermal area of Paipa is located in the eastern cordillera of Colombia, south of the Paipa town, with a surface area of 130 km2 (Figure 1). In this area, geothermal fluids rise to the surface in two specific areas; these fluids are currently only used in baths and thermal pools. Different investigations have been carried out to evaluate the geothermal potential of Paipa (Alfaro Valero C., 2002; Vásquez L., 2002; Velandia F., 2003; Cepeda & Pardo, 2004; Alfaro Valero & Espinoza, 2004; Alfaro Valero C. 2005; Alfaro Valero C., 2005b; Ortiz & Alfaro, 2010; Franco J., 2012; Vásquez L., 2012; Rodríguez & Vallejo, 2013; Llanos et al., 2015; Moyano I., 2015; Rueda J., 2016). Although several studies have assessed the geothermal system of Paipa, no research has included satellite information with thermal infrared. Thermal remote sensing has been used effectively in the exploration of geothermal resources, resulting in saving resources such as time and money, it is shown by several geothermal explorations of active volcanoes and other geothermal areas (Zhou, 1998; Ayenew, 2001; Qin et al., 2011; Norini et al., 2015; Chan et al., 2018; Darge et al., 2019). Some of these investigations have been carried out in well known geothermal areas, allowing the identification of new specific areas of geothermal interest (Coolbaugh, 2003, Gutiérrez et al., 2012). Multispectral satellites, such as Landsat 7 ETM+ and Landsat 8 OLI/TIRS, have thermal bands that are used to determine thermal anomalies in volcanic areas and geothermal areas (Mia et al., 2014, Chan & Chang, 2018, Cesarian et al., 2018; Sekertekin & Arslan, 2019). The use of thermal satellite images has some limitations; not all scenes are of a high quality due to cloudiness and spatial resolution, and due to difficult access in some areas, validation of satellite data with field data is challenging. Landsat 7 ETM+ has a thermal band with a spatial resolution of 60 m and a spectral range between 10.4 and 12.5 μm. A single-channel algorithm was applied to this band to determine surface temperature anomalies (Jiménez- Muñoz et al., 2009). Similarly, Landsat 8 OLI/TIRS has two thermal bands with a spatial resolution of 100 m and a spectral range between 10.6 to 11.19 μm (band 10) and 11.5 to 12.51 μm (band 11). A two-channel algorithm (Split Window) was applied to these two bands (Du et al., 2015). The exploration of Paipa geothermal area is being carried out for the first time by thermal satellite images. The approach used here aims to demonstrate the effectiveness of the thermal infrared satellite to determine geothermal anomalies, supported by field measurements that allow the validation of satellite data. The main objective of this research is to highlight exploratory prospects in order to the Paipa community to take advantage of new hydrothermal fluids and attain economic growth in the region. Geological and Geothermal Framework Paipa geothermal system is located in the eastern cordillera basin whose geological units are mainly sedimentary, with ages ranging from 1̴30 m.y., to the present (Figure 2). The formations of Cretaceous period correspond to Tibasosa-K1t (limestone and mudstone), Une-K1u (quartz sandstones), Churuvita group-K2ch (mudstone and sandstone), Conejo-K2c (mudstone, siltstone, and limestones), Plaeners-K2pl (siltstone and chert), Los Pinos- K2p (mudstone and siltstone), Labor-Tierna-K2lt (quartz sandstones) and Guaduas-K2E1g (mudstone, sandstones and coal). From the Paleogene period, Bogotá formation-E1-2b (quartz sandstones). From the Neogene period the Tilatá formation-N1t (sandstones and mudstones), and pyroclastic volcanic Figure 1. Location of the study area in Paipa-Boyacá, Colombia (the orange dashed rectangle marks research area). The red edge triangles represent main geothermal manifestations, and yellow edge circles represent well know geographic reference points. Finally, color polygons represent the volcanic deposits at zone. 267Determination of Hydrothermal Prospects in the Geothermal Region of Paipa (Boyacá, Colombia), Using Remote Sensing and Field Data using spectral bands of red (30m / pixel), near infrared (30m/ pixel) and thermal infrared (60m-100m/pixel) of Landsat images. We also used spectral radiance images (MOD021KM) and images of water vapor content of the product MOD05_L2 (1000m/pixel) of MODIS/ Terra, which is used to correct the thermal emissivity of the surface atmospherically. The digital elevation model- DEM (12.5m/pixel) used was obtained from the ALOS-PALSAR product, and the Pléiades image of the high-resolution visible spectrum (RGB) (0.5m/pixel) were used in conjunction with the DEM to create anaglyphs. The images cover path 7 and row 56 of Landsat images. In order to multi temporarily analyze the thermal anomalies, five Landsat diurnal images were used. These were acquired on the dates 2000-12-13, 2001-01-30, 2002-01-17, 2002-03-06, 2003-01-04 for Landsat 7 and 2016-02-01 for Landsat 8. These images were selected because they presented 0% cloudiness for the research area and are of excellent quality. Field data Like satellite images, field data is primary information. We acquired fifty field stations of solar radiation, soil, air, and surface temperatures, relative humidity, and atmospheric pressure. A radiometric station RS (powered by solar energy) and a FLIR CAT-S60 thermal camera were used. The RS (Figure 3) is composed of a data logger, a solar panel, a radiometer, an infrared radiometer, a thermistor, and a multiparameter environmental sensor. Table 1 provides details of the RS sensors. Solar radiation, air, and surface temperature data were used to correct and validate the temperature values of the satellite images. The soil temperature data was used to create a map of soil temperature anomalies in the area. Methodology This study was carried out in five stages (Figure 4). First, the reference information was established; second, the satellite images were processed deposits and intrusive domes. Also, from the Quaternary, there are colluvial (Qc), and fluvial-lacustrine (Qal) deposits located in the flat parts of the zone (conglomerates and sandstones) (Velandia, 2003). As for the volcanic deposits, these are located seven kilometers south of Paipa town (Figure 2) and are classified mainly as ash flows, ash falls, and dome blocks; the body domes composition is rhyolitic (Cepeda & Pardo, 2004). Outcropping on the surface, there are 18 geothermal manifestations, mainly hot springs and a steam vent (Figure 2) (Alfaro et al., 2017). Previous reports such as that of the Latin American Energy Organization (OLADE, ICEL, CONTECOL & Colombian Geotermia, 1982) and the Japan Consulting Institute (1983), conceptualized that this area has a geothermal interest classified from medium to high, with a possible geothermal reservoir suitable for its exploitation. Rodríguez (1998) determined that the volume of thermal water available in the area is approximately 480,000 m3. The Colombian Geological Service (CGS) proposes a source of magmatic heat, associated with the heat remaining from igneous intrusions and/or heat from radiogenic minerals, which heats the waters of the system (Alfaro Valero et al., 2017). Also, the CGS proposes the existence of two hydrothermal systems, one primary and the other secondary, for its possible temperatures in depth; the primary system would be generated by infiltration of surface water at a depth between 4 and 5 km, which would circulate from the SE of the area to the NW, reaching temperatures higher than 230°C. The secondary hydrothermal system located NW of the area, where surface water would infiltrate to a depth of 1.5 km, reaching a temperature around 75°C, restricted laterally by the Canocas and Hornitos faults, superficially this system does not show hot springs. Materials and methods Satellite Data Satellite images are primary information sources in this investigation. These images include scenes from Landsat 7, Landsat 8, MODIS/Terra, ALOS- PALSAR and Pléiades. Land surface temperature anomalies are determined Figure 2. Geologic Map of the study area, from Velandia (2003) and Cepeda & Pardo (2004). 268 Rafael Andrés Calderón-Chaparro, German Vargas-Cuervo digitally; third, index maps were defined; fourth, a spatial model was applied to the four index maps, and finally, the hydrothermal prospects were defined. Digital image processing This stage consists of three steps. First, the Landsat bands were radiometrically calibrated and atmospherically corrected. Second, the red and infrared near bands were used to calculate the thermal emissivity of the surface. Finally, the algorithms of a single thermal band (Jiménez-Muñoz & Sobrino, 2003) and two Split-Window thermal bands (Du et al., 2015) were used to determine the thermal anomalies of the surface. The following sections provide a detailed description of each of these procedures. Land Surface Temperature retrieval The first step to obtaining surface temperature images is the conversion of digital numbers (DN), in spectral radiance values (L𝜆) to the top of the TOA atmosphere (Chander et al., 2009). For Landsat 7 and Landsat 8 images equations (1) USGS (2018) and (2) USGS (2016) were used. L G Q Brescale Cal rescale = +*� � (1) L M Q AL cal L = +* � (2) Where, Grescale, Brescale, ML and AL, are specific rescaling factors for each band; they are found in the metadata of each of the images. QCal is the raw value of each pixel in the images, for the case of Landsat 7, those values are between 0 and 255, and for the Landsat 8 images values are between 0 and 65,535. Atmospheric Correction Atmospheric effects correct the bands in the visible and near-infrared spectrum. This correction is applied due to absorption and dispersion processes, as a result of the existence of water vapor, gases (Ozone, Oxygen, Carbon Dioxide, Nitrogen Oxide, etcetera.) and molecules in the atmosphere, which affect the amount of radiation received by the sensor. These atmospheric effects are removed in order to obtain an image of spectral reflectance whose values reflect real physical magnitudes of the surface. To obtain atmospherically corrected images, the atmospheric correction module Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) was used. In general, FLAASH is an atmospheric correction method based on the MODTRAN 4.0 Radiative Transfer code, which allows the removal of the effects of dispersion and absorption in the Multispectral and Hyperspectral images, which was developed by the Laboratory of the US Air Force. (Kayadibi, 2011). Thermal emissivity calculation Once the red and near-infrared bands were corrected atmospherically (Example: bands 3 and 4 for Landsat 7), the values of surface emissivity were determined. The approach used here is the one proposed by Sobrino et al., (2004), it is based on the determination of the emissivity for each pixel, based on the proportion of vegetation Pv, obtained from the values of the Normal Differential Vegetation Index (NDVI), equation (3). NDVI R R R R nir Red nir Red = − + (3) The method is called NDVI thresholds (NDVITH), which considers three different surface classes, dependent on the NDVI values: a) Bare soil (NDVI <0.2), b) mixed surface (0.2≥NDVI≤0.5), and c) Surface fully covered with vegetation (NDVI> 0.5). a) Bare soil (NDVI <0.2): The emissivity value for Landsat 7 images is established as constant ε = 0.97. In the case of Landsat 8 images, the emissivity was determined from the reflectance of the red band ρred, for each of the thermal bands, ε ρ10 0 973 0 047= −. . red and ε ρ11 0 984 0 026= −. . red (Yu et al., 2014). b) Mixed surface (0.2 ≥ NDVI ≤ 0.5): The pixels are considered as a mixture between vegetation and bare soil. The emissivity is obtained as a function of the proportion of vegetation, as follows,    = + −( ) +v v s vP P d1 , where, εv is the emissivity of the vegetation, εs is the soil emissivity, and dε is a term which includes the effect of geometric distribution of natural surfaces and internal reflections (Sobrino et al., 2004). The proportion of vegetation is calculated from the following equation, P NDVI NDVI NDVI NDVIv min max min = − −       2 , where, NDVImin = 0 2. and NDVImax = 0 5. . dε is calculated in the following way d P Fs v v  = −( ) −( )1 1 * * , where F is a form factor, with a mean value of 0.55, which considers different kinds of geometric distribution. Figure 3. Radiometric Station (RS) used in field acquisition. Table 1. Instruments used in the field data collection Variable Units Sensor Brand / Model Soil Temperature °C Thermistor Apogee / ST100 Surficial Temperature °C Non-contact Infrared thermometer Apogee / SI-4H1 Object surficial Temperature °C Thermal camera CAT FLIR / S60 Air Temperature °C VP-4 Decagon / PASSRHT Vapor pressure kPa Relative Humidity % Atmospheric Pressure kPa Incident and reflected solar radiation W/m2 Net Radiometer Apogee / SN500 Net Radiation W/m2 269Determination of Hydrothermal Prospects in the Geothermal Region of Paipa (Boyacá, Colombia), Using Remote Sensing and Field Data For the case of the Landsat 8 thermal bands, the average emissivity values of the soil and vegetation of each one were taken into account: εs10= 0.9668, εs11= 0.9747, εv10= 0.9863, εv11= 0.9896 (Yu et al., 2014). c) Surface fully covered with vegetation (NDVI> 0.5): It was assumed for both Landsat 7 and Landsat 8 images, a constant value ε = 0.99. Calculation of Surface Kinetic Temperature Two algorithms were chosen to determine the kinetic surface temperature of Paipa’s geothermal area. For the Landsat 7 ETM + thermal band, the one- channel algorithm proposed by Jiménez-Muñoz & Sobrino (2003) was used, and for the Landsat 8 TIRS thermal bands, the Split-Window algorithm, proposed by Du et al., (2015) was used. These two algorithms are based on the radiative transfer equation (RTE), equation (4). This equation expresses the law of Planck, defining that the spectral radiance B(λ,T) of a black body is a function of the wavelength λ(m) of the electromagnetic energy and the absolute temperature of the body T(K). B T c exp c T    ,( ) =     − 1 5 2 1 (4) Where c1 and c2 are constants of Planck, with values of 1.19104 * 108 W.µm4.m-2.sr-1 and 14387.7 μm.K, respectively. Single Channel Algorithm (SC) Jiménez-Muñoz & Sobrino (2003) and Jiménez-Muñoz et al., (2009) proposed an algorithm to determine the surface temperature LST, from a single channel, equation (5), based on the concept of Atmospheric Functions (FAts), which depend on the atmospheric transmissivity and the descending and ascending atmospheric radiations. LST Lsen= +( ) +    +γ ε ϕ ϕ ϕ δ 1 1 2 3 (5) Where ε is the emissivity of the surface; γ and δ are two parameters dependent on the Planck functions, expressed in the following way,   ≈ T b L sen sen 2 and δ γ ≈ −T T bsen sen 2 , Where Tsen is the temperature of brightness in the sensor (equation 6); bγ is a constant, with a value of 1277 K for Landsat 7 and Lsen is the spectral radiance in the sensor. Finally, φ1, φ2, φ3 are atmospheric functions, ϕ ω ω1 20 14714 0 15583 1 1234= − +. . . ; ϕ ω ω2 21 1836 0 37607 0 52894= − − −. . . ; ϕ ω ω3 20 04554 1 8719 0 39071= − + −. . . ϕ ω ω3 20 04554 1 8719 0 39071= − + −. . . , where ω is the water vapor content (WVC) in the atmosphere in g/ cm2. The product MOD05_L2 of the MODIS sensor, offers continuous values of ω at a spatial resolution of 1 km2, being appropriate due to the small spatial variation of ω in the atmosphere (Ren et al., 2015). T K ln K L sen = +     2 1 1  (6) Where, K1 [W.m-2.sr-1.µm-1] and K2 [K] are calibration constants (Table 2); Lλ is the spectral radiance in the sensor [W/(m2.sr.µm)] and ln is the natural logarithm. Before being able to use the MODIS images, a process of georeferencing of the MOD021KM products was carried out, starting from the Landsat images. Once the MOD021KM images were georeferenced, a co-registration process was performed with their respective MOD05_L2 images, obtaining correctly geo-positioned WVC images. Once the WVC images were correctly geo- positioned, their values were validated using the values of relative humidity and air temperature, collected in the field, by applying equation 7, where T0 is Figure 4. General Methodological Framework. 270 Rafael Andrés Calderón-Chaparro, German Vargas-Cuervo the surface air temperature, RH is the relative humidity of the air and WVC is given in g/cm2. WVC exp T T RH= +               0 0981 10 0 6108 17 27 237 3 0 0 . * * . * . * . *  + 0 1697. (7) The WVC values obtained from the field data vary from 0.52 g/cm2 to 1.46 g/cm2, which are very similar to the average values of WVC obtained by the MODIS sensor in each of the dates of the Landsat images (Table 3). This similarity allowed the validation of the WVC satellite values. Table 3. Water vapor content statistics of the MOD05_L2 images, for each of the dates used in Landsat. Statistics MOD05 20001213 MOD05 20010130 MOD05 20020117 MOD05 20020306 MOD05 20030104 Min 0.577282 0.486306 0.369391 0.706861 0.955387 Max 1.54891 0.93637 1.032812 1.71515 2.12305 Mean 1.147071 0.651245 0.749141 1.145315 1.33665 Split Window Algorithm (SW) Du et al., (2015) proposed a practical algorithm, which takes advantage of the two thermal bands of the TIRS sensor, onboard the Landsat 8 satellite. The algorithm obtains the surface temperature from a linear relationship, with the combination of the brightness temperatures of the two thermal bands, with clear sky conditions. Like the algorithm of Jiménez-Muñoz & Sobrino (2003), the SW is based on the theory of radiative transfer, so it is necessary to have atmospheric information at the time of image acquisition. To solve this, the algorithm of equation (8) depends on eight coefficients, which vary depending on the WVC, at the moment of image acquisition. The coefficients were calculated from a modified SW method of covariance-variance relation, obtaining five different ranges of WVC, [0.0 - 2.5], [2.0 - 3.5], [3.0 - 4.5], [4.0 - 5.5], and [5.0 - 6.3] g/cm2. LST b b b b T T b b bi j= + + − +     + + + − +    0 1 2 3 2 4 5 6 2 1 2 1        ∆ ∆ TT T b T Ti j i j + + −( )2 7 (8) Where Ti and Tj are the brightness temperatures of bands 10 and 11, respectively; ε is the average emissivity of the thermal bands    = +i j 2 ; Δε is the difference of the emissivity Δε=εi-εj, and b0, b1, ...b7 are the coefficients determined for each WVC range. The coefficients used correspond to the values defined for a range of WVC between 0.0 and 2.5 g/cm2, determined thanks to the record of the MODIS images and the field values, whose average value is 1.3 g/cm2. Obtaining Structural Lineaments The first step to obtain the structural lineaments, was to combine the ALOS-PALSAR (DEM) scenes in a single image; for this, the mosaic tool of ERDAS Imagine 2018® was used. Subsequently, shadow models were obtained; these models were combined with the DEM and the Pléiades images of the visible spectrum, thus obtaining three-dimensional images (anaglyphs). Once the anaglyphs were obtained, the structural analysis of the geothermal area of Paipa was carried out. The deformation structures were interpreted visually, using the 3D views of the anaglyphs. The lineaments were mapped using a geographic information system (GIS). Moreover, the identified structures were verified with the existing geological information (Velandia, 2003) and the fieldwork carried out in this study. Iso-resistivity map Based on the reference information, the iso-resistivities map published by the CGS (Franco, 2016 in Alfaro Valero et al., 2017) was used. This iso- resistivity surface information was measured at an elevation of 2,500 m.a.s.l., and the anomalies of low resistivity indicate the location of geothermal fluids in the area. Soil Temperature Map A soil temperature map of the Paipa geothermal system was made. One hundred ninety-seven temperature probes were used, recorded at depths greater than one meter. The Geostatistical wizard tool of ArcMap© was used by applying an interpolation method of kriging. The one hundred ninety-seven soundings are constituted by 147 soundings taken by the CGS (Rodríguez & Vallejo, 2013) and 50 new temperature data acquired in the field throughout this study. Spatial Model The maps obtained from satellite surface temperature (Sst), structural lineaments (Sl), iso-resistivity surface (Is), and soil temperature (St) were integrated in an unweighted way into a spatial model. Three index values were defined to each of these maps, corresponding to the areas of higher temperature (SST and St), areas of influence of the structural lineaments (Sg) and areas of lower resistivity (Is). The model was based on obtaining a spatial relationship between the four index maps. To achieve this, the Union spatial function, offered by ArcMap©, was used. The result obtained called Integrated Index Map (Ii), corresponds to new polygons, containing in their attributes the values of the indexes that were superimposed (Figure 5). Finally, the index values in each polygon were summed and divided into the total number of indexes (Equation 8). Ii I I I I n Tss Le Ir Ts= + + +� � � (8) Figure 5. Illustrative diagram of the result obtained from the Union process. The input data is three layers, each with a defined attribute value (1, 2 and 3); the result presents the same figure, but now there are seven different polygons, and each contains a different number of attributes; note the orange polygon now contains the attribute values 1, 2, and 3. Index Maps The index maps correspond to geophysical and geological indicators, which have a spatial relationship with the occurrence of hydrothermal fluids on the surface; each index was classified into three qualitative ranges depending on their temperature, resistivity and lateral area of influence of the lineaments. The Table 2. Calibrations constants used to convert spectral radiance in brightness temperature, for Landsat 7 ETM+ and Landsat 8 TIRS. Satellite Band K1 K2 Landsat 7 6 666.09 1282.71 Landsat 8 10 774.8853 1321.0789 11 480.8883 1201.1442 271Determination of Hydrothermal Prospects in the Geothermal Region of Paipa (Boyacá, Colombia), Using Remote Sensing and Field Data defined index values are 1, 2, and 3, classified qualitatively as, low, medium, and high probability, respectively. a) Surface temperature index map Two of the thermal satellite images were integrated. The images used correspond to a Landsat 7 ETM+ of March 6th, 2002 and a Landsat 8 TIRS of February 1st, 2016. These images were selected for the dominant surface coverage on those dates. The ground covers correspond to surfaces mainly of half-bare soil, and this is ideal since the temperature values can be affected by the evapotranspiration processes of the vegetation. Thus, the resulting temperature anomalies of the images are similar in shape and magnitude. Table 4 presents the index values defined for the three highest temperature ranges of satellite thermal images. Table 4. Satellite temperature indices. Three ranges of maximum temperature were defined and classified qualitatively. Image Index Temperature Range (°C) Qualitative Classification Landsat 7 y 8 1 39 – 41 Low 2 41 – 43 Medium 3 43 - 50 High The temperature anomalies of the two Landsat images were classified according to Table 4, later, they were converted into a vector format, and the ArcMap© union process was applied. The result obtained corresponds to the integrated temperature anomalies of the two satellite images (Figure 12-a). b) Lineaments Index Map Based on the mapped lineaments, three influence areas (buffer) were established around them. These three areas correspond to the three established rates, low, medium, and high probability of allowing the circulation of hydrothermal fluids (Table 5). The index map of structural lineaments is presented in Figure 12-b. Table 5. Indices of lineaments. Three influence distances were defined with their respective qualitative classification. Index Influence Range (m) Qualitative Classification 1 50 Low 2 100 Medium 3 150 High c) Iso-resistivity index map From the iso-resistivities surface (figure 10), the index values of this map were defined. The values of resistivity found in the zone vary between 0 and 2,251 Ω.m, and according to the geothermal conditions of Paipa, it was interpreted that the low values of resistivity between 0 and 17.3 Ω.m, are the result of hot saline fluids (Ussher et al., 2000). These low values of resistivity were taken and classified in the three index ranges (Table 6), defining the iso-resistivity index map (Figure 12-c). Table 6. Indices defined for three minimum ranges of subsoil resistivity, with their respective qualitative classification. Index Resistivity Range (Ohm.m) Qualitative Classification 1 11.6 – 17.3 Low 2 7.7 – 11.6 Medium 3 < 7.7 High d) Soil temperature index map A frequency analysis of the soil temperature data was carried out, concluding that the majority of the data set corresponds to values without geothermal influence, that is, the values of temperature between 15 °C and 19 °C. On the other hand, it was established that the temperature data between 19.1 °C and 37 °C correspond to values with geothermal influence (Figure 6). From this last temperature range, the index values were established (table 7), and the soil temperature index map was obtained (Figure 12-d). Prospects The prospects of the Paipa geothermal area were obtained by making a spatial superposition of the four index maps. The objective was to identify the areas where the highest index values were superimposed (index value 3), that means, high satellite temperature, high soil temperature, low electrical resistivity and proximity to the fault zone. Results and discussion Land Surface Temperature Anomalies The results of the Landsat bands processing are presented in figure 7. The application of the algorithms proposed by Jiménez-Muñoz & Sobrino (2003) and Du et al., (2015) in the thermal bands of Landsat 7 ETM+ and Landsat 8 TIRS, allowed to correct in these, the effect of the emissivity of the surfaces and the processes of dispersion and absorption generated by atmospheric conditions, mainly due to the effect of the water vapor content in the atmosphere. The correction of WVC in the atmosphere was carried out satisfactorily with the MOD05_L2 products, for the case of the “Single Channel” algorithm and the case of the “Split Window” algorithm with the chosen range of water vapor content between 0.0 and 2.5 g/cm2. This was due to the atmospheric data of temperature and relative humidity of the air, acquired in the field, which allowed it to validate the research areas with an average WVC value of 1.3 g/ cm2, validating the values of WVC recorded by the MODIS sensor and the chosen range for the “Split window” algorithm. Another relevant aspect of the application of the algorithms is obtaining surface temperature anomalies with higher temperature magnitudes, and the reduction in the size of the anomalies in comparison with the temperature anomalies obtained without the atmospheric and emissivity corrections (Temperature of brightness). As an example, Figures 8a and 8b are compared, and it is possible to observe the differences in the size of the temperature anomalies, being of much larger extension for the temperature anomalies of brightness, while the temperature anomalies where the corrections are applied have a smaller size and their magnitude (temperature) is higher. Regarding the variation of the temperature magnitudes, there was a generalized increase Fr eq ue nc y Temperature Range ºC 80 70 60 50 40 30 20 10 0 [1 5, 1 6] [1 6, 1 7] [1 7, 1 8] [1 8, 1 9] [1 9, 2 0] [2 0, 2 1] [2 1, 2 2] [2 2, 2 3] [2 3, 2 4] [2 4, 2 5] [2 5, 2 6] [2 6, 2 7] [2 7, 2 8] [2 8, 2 9] [2 9, 3 0] [3 0, 3 1] [3 1, 3 2] [3 2, 3 3] [3 3, 3 4] [3 4, 3 5] [3 5, 3 6] [3 6, 3 7] Figure 6. Histogram of frequencies of soil temperatures in the Paipa area. Table 7. Indices defined for three maximum ranges of soil temperature, with their respective qualitative classification. Index Temperature Range (°C) Qualitative Classification 1 19 – 22 Low 2 22 – 25 Medium 3 25 – 37 High 272 Rafael Andrés Calderón-Chaparro, German Vargas-Cuervo Figure 7. Anomalies of land Surface temperature from Landsat 7 ETM+ on December 13, 2000; January 30, 2001; January 17, 2002; March 06, 2002; January 04, 2003; and Landsat 8 TIRS on February 01, 2016. 273Determination of Hydrothermal Prospects in the Geothermal Region of Paipa (Boyacá, Colombia), Using Remote Sensing and Field Data Figure 8. Comparative images between the satellite thermal anomaly sizes, from the Landsat 7 EMT+ scene, of March 6, 2002. a) Brightness temperature image. b) Surface temperature image. Figure 9. Structural lineaments map of Paipa geothermal zone. 274 Rafael Andrés Calderón-Chaparro, German Vargas-Cuervo of around +8°C, for the hottest anomalies, while for the colder temperature anomalies there was an increase of around 3.5°C. The satellite thermal anomalies reflect an excess of heat coming from the interior of the earth. Thus, in areas of geothermal activity, the images of surface temperature that are obtained have two sources of energy, one coming from solar heating and another coming from the interior of the earth. The equation of the surface energy balance (equation 9) allows us to comprehend how it is possible to identify zones with geothermal activity. Q H E Gn = + +�  (9) Where Qn is the net solar radiation received by the surface; H is the sensible heat flow between the ground floor interface; λE is the latent heat flux in the transition phase of water between the surface and the atmosphere; G is the heat flow of the soil. Locally, it is assumed that H and λE have little influence on surface thermal anomalies so that the net radiation and the heat flow from the soil are those that would generate the signature of heat registered by the sensor. Structural Lineaments Map Figure 9 presents the map of structural lineaments obtained from the 3D interpretation. The result shows the existence of a dense network of fractures which would affect the geothermal area of Paipa. The lineaments validated in the field and reference information corresponds to the Hornito, Canocas, El Batán, El Bizcocho, Cerro Plateado, and Las Peñas faults. Iso-resistivity Map Figure 10 shows the iso-resistivities surface of the Paipa geothermal area. Resistivities under 17.3 Ω.m were interpreted as a result of a geothermal source. Six specific zones of high conductivity (<2.3 Ω.m) were identified in the study area. The anomaly 1 in the north of the area, is not taken into account because this anomaly is generated by a deposit of salt in that area. The other anomalies are interpreted as shallow accumulations of geothermal fluids. Soil Temperature Map Figure 11 shows the soil temperature map, with data acquired at one- meter depth. The presence of six anomalous zones is clearly identified, with two principal anomalies, one located to the east of Sochagota lake (anomaly 1) and the second to the south, in the sector of the El Delfín and La Playa pools (anomaly 4), each corresponding spatially to the zones where the geothermal manifestations of the zone are presented. In general, there is a spatial correlation between resistivity anomalies and soil temperature anomalies, mainly between anomalies number 2 and number 1 (Figure 10 and 11), anomalies 3-4-5 with anomaly 4 (Figure 10 and 11), and the number 6 of the two anomalies (Figure 10 and 11). Spatial Model - Index maps The implemented spatial model took into account indicators with a spatial relationship with the presence of hydrothermal fluids (Figure 12). Lineaments act as migration paths or may restrict the lateral distribution of hot fluids. The resistivity anomalies reflect the possible zones of accumulation of geothermal fluids because these fluids have physical and chemical properties characteristic of temperature and salinity, which produce very high conductivities. Soil temperature anomalies at shallow depths are the direct manifestations of heat transfer from the subsoil to the surface, with convection being its transfer mechanism. Finally, the thermal anomalies obtained from the satellite images are the continuous representation of the surface temperatures mainly product of the interaction between the solar radiation and the heat coming from the interior of the earth, which allows identifying anomalously hot zones through the surface. These variables integrated from a geographic information system, in an unweighted way, facilitated the spatial correlation between each of them and allowed for the application of a simple algorithm with which it was possible to generate point zones, which are assumed as the areas of most significant probability to find new hydrothermal sources (prospects). These results allow us to affirm that the defined methodology is reliable, meets the stated objective, and can be applied to any geothermal zone, where the same indicators used for the model are available. Prospects The areas of exploratory interest in the geothermal area of Paipa are presented in Figure 13-b. The results showed six areas with hydrothermal potential, of which four were defined as exploratory prospects, and two correspond to hydrothermal springs currently in production (hot springs of the ITP Paipa Tourism Institute, and El Delfín and La Playa pools). These results validated the applied methodology and can be used in any geothermal zone, where the same geological-geophysical indicators used for this model are available (satellite temperature maps, structural lineaments, geoelectric studies, and soil temperature data). The four prospects were listed from one to four, and are located in: 1) El Durazno; 2) East of Alto Los Volcanes, divided into two circular anomalies; 3) Alluvial valley of the Honda creek, and 4) East of the El Delfín-La Playa pools. Validation of Results Surface Temperature Surface temperature anomalies can be used with confidence once they have been appropriately validated. The approach used here includes the spatial comparison between surface temperature anomalies and soil temperature anomalies since the latter are a direct expression of the heat present in the subsoil. In general, five of the six soil temperature anomalies correspond spatially with satellite surface temperature anomalies (Figure 14 a-1), demonstrating the geothermal component in thermal satellite images. Exploratory Prospects The results showed six areas of exploratory interest, two of which coincide spatially with the superficial geothermal manifestations of Paipa (Figure 13b). This spatial coincidence validates the results of the implemented method. The production area #1 (figure 15), is located in the main discharge area of the system, which contains approximately ten hot springs, with temperatures between 40°C and 72.1°C. These thermal waters are currently being used by the center Tourism (Instituto de Turismo de Paipa, ITP) and hotels in the area. Production area #2 (figure 16) corresponds with the second discharge area of the system. In this area, there are eight geothermal manifestations, of which seven are hot water springs and one steam vent (Alfaro Valero et al., 2017). The thermal waters register temperatures up to 60°C and are used in the El Delfín and La Playa pools. The steam vent registers a temperature of 64.6°C. Geological Context Unlike geothermal reservoirs in active volcanic zones where thermal anomalies are located following a linear pattern according to the fault lineaments (Wu et al., 2012, Norini et al., 2015), the geothermal system of Paipa is different as the geothermal reservoirs are located through sedimentary rocks, which generates a lateral distribution of the hot fluids. The location of the thermal anomalies is found mainly to the northwest of the geothermal area (Figure 17a). This location can be a product of the circulation of hot fluids at shallow depth, through porous stratigraphic units of the Tilatá, Bogotá, Guaduas, Labor- Tierna, and Los Pinos formations, and facilitated by the faults, which function as migration paths. Figures 17a and 17b show the El Batán overthrust fault, 275Determination of Hydrothermal Prospects in the Geothermal Region of Paipa (Boyacá, Colombia), Using Remote Sensing and Field Data Figure 10. Iso-resistivity anomalies map of Paipa geothermal zone. From Franco (2016) in Alfaro et al. (2017). Figure 11. Soil temperature map of Paipa geothermal zone. 276 Rafael Andrés Calderón-Chaparro, German Vargas-Cuervo Figure 12. Index maps of Paipa geothermal zone. a) Land surface temperature index-map, b) Structural Lineaments index map, c) Iso-resistivity Index Map, d) Soil temperature Index Map. 277Determination of Hydrothermal Prospects in the Geothermal Region of Paipa (Boyacá, Colombia), Using Remote Sensing and Field Data Figure 13. a) Integrated Index Map of Paipa Geothermal Zone. b) Explorations Prospects of Paipa Geothermal zone. with a northeast course, of vergence towards the east, functioning as a limit (restricting), for the lateral circulation of hot fluids. The anomalies of surface temperatures found between the El Hornito and Canocas faults, southwest of Lake Sochagota (Figure 17a), would be validated, according to the conceptual geothermal model proposed by the CGS (Alfaro Valero et al., 2017), who proposed the existence of a secondary thermal system differing from the main one. This secondary system has a circulation of hot fluids at a shallow depth towards the northeast, by the porous sediments of the Guaduas and Bogotá formations. Differently, for the thermal anomalies south of the Canocas fault, between the corridor generated by the faults the Bizcocho and the Batan, can be validated by the presence of hydrothermal fluids, inferred by the low resistivity anomalies (<1 Ω.m), obtained from the magnetotelluric technique at a height of 2,450 m.a.s.l, see figure 17b (Alfaro Valero et al., 2017). Other satellite thermal anomalies of smaller size (some coinciding with the Cerro Plateado fault on the alluvial deposits and others on volcanic deposits, east of the Alto Los Godos, coinciding with the Godos Lineament) have some spatial coincidence with low resistivity anomalies of the Magnetotelluric technique (Figure 17b). Finally, the satellite thermal anomalies present in the ITP-Infiboy and El Delfín-La Playa sectors are validated because they are present in the two zones where the only surface geothermal manifestations of the system are presented. These manifestations are produced through alluvial deposits, for the ITP-Infiboy sector facilitated by the crossing of the El Hornito and El Bizcocho faults, and for the El Delfín-La Playa sector, facilitated by the El Batán fault (Alfaro Valero et al., 2017). Conclusions This research used thermal infrared technology combined with conventional exploration methods to identify hydrothermal prospects in the Paipa geothermal region. Surface temperatures were obtained from the research area, from Landsat 7 ETM+ and Landsat 8 OLI/TIRS images. These images were calibrated and atmospherically corrected using the single channel and Split window algorithms, which led to spatially improved temperature anomalies being obtained. Additionally, fieldwork was carried out to validate satellite thermal anomalies. The fact that most satellite thermal anomalies coincide spatially with soil temperature anomalies validates the utility of satellite thermal bands in geothermal exploration. The spatial model implemented, determined six prospective areas of which two were validated by the spatial correlation of surface geothermal manifestations of the geothermal area of Paipa. The use of thermal satellite images for the correct identification of geothermal anomalies requires the incorporation of geological-geophysical analysis; this integration considerably improves the interpretation of satellite thermal information. 278 Rafael Andrés Calderón-Chaparro, German Vargas-Cuervo Figure 14. Spatial comparison of soil temperature anomalies and satellite temperature anomalies. a-b) Anomalies to the southeast of Lake Sochagota; c-d) Anomalies to the southwest of Lake Sochagota towards the Durazno sector; e-f) Anomalies to the sector of El Delfín and La Playa pools, where the presence of the second most important manifestations of the Paipa geothermal system occurs. 279Determination of Hydrothermal Prospects in the Geothermal Region of Paipa (Boyacá, Colombia), Using Remote Sensing and Field Data Continuation figure 14. Spatial comparison of soil temperature anomalies and satellite temperature anomalies. g-h) Anomalies to the southeast of the El Delfín and La Playa pools, on the alluvial valley of the Honda creek; i-j) Anomalies to the east of the El Delfín and La Playa pools, close to Cemex quarry; k-l) Anomalies to the southwest of the El Delfín and La Playa pools, and northwest of Los Volcanes high. 280 Rafael Andrés Calderón-Chaparro, German Vargas-Cuervo Figure 15. a) The spatial location of prospect #1 and hot springs in the main discharge area of the system (Infiboy-ITP). b) Thermal and visible photography of the Hot spring located next to the Hotel Lanceros; a maximum temperature of 51.3°C was recorded. c) Thermal and visible picture from the hot spring of the devil’s eye; a temperature of 72.1°C was registered. Figure 16. a) The spatial location of prospect #3 and hot springs in the second discharge area of the system (El Delfin-La Playa). b) Thermal and visible picture of the spring III El Delfín, it has a temperature of 54.1°C. c) Thermal and visible picture of hot springs I and II El Delfín. d) Thermal and visible picture of the steam vent in the La Playa sector; a maximum temperature of 64.6°C was registered. 281Determination of Hydrothermal Prospects in the Geothermal Region of Paipa (Boyacá, Colombia), Using Remote Sensing and Field Data Figure 17. a) The spatial relationship between satellite thermal anomalies of Landsat 7 ETM+, from March 6, 2006, and surface geological formations; b) Spatial relationship between the anomalies of low resistivity (Magnetotelluric, Alfaro et al., 2017) and the lineaments of the research area. Acknowledgments This investigation was developed within the framework of the 733 calls, from the Government of Boyacá and Colciencias. We thank the Government of Boyacá for sponsoring the entire project, as well as the National University of Colombia for its administrative collaboration in the management of resources. We also thank NASA and the Alaska Satellite Facility for providing the Landsat images and digital elevation models used in this investigation. References Alfaro Valero, C. M. (2002). Estudio Isotópico de Aguas del Área Geotérmica de Paipa. Colombian Geological Survey (Previously INGEOMINAS). Subsoil Resources Branch. Bogotá, Colombia. Alfaro Valero, C. M. (2005). Alteración hidrotermal en el sistema geotérmico de Paipa. Colombian Geological Survey (Previously INGEOMINAS). Subsoil Resources Branch. Bogotá, Colombia. Alfaro Valero, C. M. (2005). Geoquímica preliminar de gases del sistema geotérmico de Paipa. Colombian Geological Survey (Previously INGEOMINAS). Subsoil Resources Branch. Bogotá, Colombia. Alfaro Valero, C. M., & Espinoza, O. (2004). Sondeo preliminar de radón en el área geotérmica de Paipa. Colombian Geological Survey (Previously INGEOMINAS). Subsoil Resources Branch. Bogotá, Colombia. Alfaro Valero, C. M., Matiz León, J. C., Rueda Gutiérrez, J. B., Rodríguez Rodríguez, G. F., Gonzáles Idárraga, C. E., Beltrán Luque, M. Á., . . . Malo Lázaro, J. E. (2017). Actualización del Modelo Conceptual del Área Geotérmica de Paipa. Colombian Geological Survey. Basic Geosciences. Bogotá, Colombia. Ayenew, T. (2001). Surface Kinetic Temperature Mapping Using Satellite Spectral Data in Central Main Ethiopian Rift and Adjacent Highlands. Ethiopian Journal of Science, 24, 51-68. Cepeda, H., & Pardo, N. (2004). Vulcanismo de Paipa. Colombian Geological Survey (Previously INGEOMINAS). Subsoil Resources Branch. Bogotá, Colombia. Cesarian, D. P., Abir, I. A., & Isa, M. (2018). Comparison of In-Situ Temperature and Satellite Retrieved Temperature in Determining Geothermal Potential in Jaboi Field, Sabang. Journal of Physics: Conference Series. Indonesia, December, 032017. Chan, H.-P., & Chang, C.-P. (2018). Exploring and monitoring geothermal and volcanic activity using Satellite Thermal Infrared data in TVG, Taiwan. Terrestrial Atmospheric and Oceanic Sciences, 29 (4), 387-404. Chan, H.-P., Chang, C.-P., & Dao, P. D. (2018). Geothermal Anomaly Mapping Using Landsat ETM+ Data in Ilan Plain, Northeastern Taiwan. Pure and Applied Geophysics, 175, 303-323. Chander, G., Markham, B. L., & Helder, D. L. (2009). Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors. Remote Sensing of Environment, 113, 893-903. Coolbaugh, M. (2003). The Prediction and Detection of Geothermal Systems at Regional and Local Scales in Nevada using a Geographic Information System, Spatial Statistics, and Thermal Infrared Imagery. Ph.D. Thesis, Department of Geological Sciences and Engineering, University of Nevada, Reno, United States of America. Corpoboyacá. (2006). Plan de ordenación y manejo ambiental de la cuenca alta del río Chicamocha. Technical report, Pedagogical and Technological University of Colombia, Tunja, Boyacá. Darge, Y. M., Hailu, B. T., Muluneh, A. A., & Kidane, T. (2019). Detection of geothermal anomalies using Landsat 8 TIRS data in Tulu Moye geothermal prospect, Main Ethiopian Rift. International Journal Applied Earth Observation and Geoinformation, 74, 16-26. Dickson, M., & Fanelli, M. (2003). Geothermal Energy, Utilization, and Technology. UNESCO, Paris, France, 221 pp. DiPippo, R. (2015). Geothermal Power Plants: Principles, Applications, Case Studies, and Environmental Impact (Fourth ed.). Butterworth- Heinemann print, North Dartmouth, Massachusetts, United States of America, 800 pp. Du, C., Ren, H., Qin, Q., Meng, J., & Zhao, S. (2015). A Practical Split-Window Algorithm for Estimating Land Surface Temperature from Landsat 8 Data. Remote Sensing Journal, 7, 647-665. Franco, J. (2012). Prospección Geoeléctrica en los Alrededores del Volcán de Paipa-Boyacá. Colombian Geological Survey. Research and Exploration of Geothermal Resources, Subsoil Resources Branch, Bogotá, Colombia. Gutiérrez, F., Lemus, M., Parada, M., Benavente, O., & Aguilera, F. (2012). Contribution of ground surface altitude difference to thermal anomaly detection using satellite images: application to volcanic/geothermal complexes in the Andes of Central Chile. Journal of Volcanology and Geothermal Research, 237-238, 69-80. 282 Rafael Andrés Calderón-Chaparro, German Vargas-Cuervo Japan Consulting Institute (1983). Feasibility study report of geothermal power plant for Instituto Colombiano de Energía Eléctrica. Jiménez-Muños, J. C., & Sobrino, J. A. (2003). A Generalized Single-Channel Method for Retrieving Land Surface Temperature from Remote Sensing Data. Journal of Geophysical Research, 108, 1-9. Jiménez-Muñoz, J. C., Cristóbal, J., Sobrino, J. A., Sòria, G., Ninyerola, M., & Pons, X. (2009). Revision of the Single-Channel Algorithm for Land Surface Temperature Retrieval from Landsat Thermal-Infrared Data. IEEE Transactions on Geoscience and Remote Sensing, 47(1), 339-349. Kayadibi, Ö. (2011). Evaluation of Imaging Spectroscopy and Atmospheric Correction of Multispectral Images (ASTER and Landsat 7 ETM+); 5th International Conference on Recent Advances in Space Technologies - RAST2011, IEEE Xplore, Istanbul, Turkey, July, 154-159. Llanos Rodríguez, E. M., Bonet, C., & Zengerer, M. (2015). 3D Geological - Geophysical Model Building and Forward and Inverse Modeling of Magnetism and Gravimetry Data from Paipa Geothermal Area, Colombia - Final Report. Melbourne. Lund, J., & Boyd, T. (2015). Direct utilization of Geothermal Energy 2015 Worldwide Review, Proceedings World Geothermal Congress, Melbourne, Australia, April, 1-31. Mia, M. B., Nishijima, J., & Fujimitsu, Y. (2014). Exploration and monitoring geothermal activity using Landsat ETM + images A case study at Aso volcanic area in Japan. Journal of Volcanology and Geothermal Research, 275, 14-21. Moyano, I. (2015). Modelamiento de Estructuras Geotérmicas Ocultas a partir de la Interpretación de Anomalías Magnetotelúricas de Resistividad en la Región de Paipa, Boyacá-Colombia. M.Sc. Thesis, Department of Geosciences, National University of Colombia, Bogotá, Colombia. Norini, G., Groppelli, G., Sulpizio, R., Carrasco-Núñez, G., Dávila-Harris, P., Pellicioli, C., . . . De Franco, R. (2015). Structural analysis and thermal remote sensing of the Los Humeros Volcanic Complex: Implications for volcano structure and geothermal exploration. Journal of Volcanology and Geothermal Research, 301, 221-237. OLADE, ICEL, CONTECOL, & Geotermia Colombiana. (1982). Report, Estudio de reconocimiento de los recursos geotérmicos de la república de Colombia. Ortiz, I. D., & Alfaro, C. (2010). Inventario de Puntos de Agua y Geoquímica de las Áreas Geotérmicas de Paipa e Iza: Aguas, Suelos y Peloides. Colombian Geological Survey. Exploration of Geothermal Resources, Subsoil Resources Branch, Bogotá, Colombia. Qin, Z., Zhang, N., Nan, P., & Chai, L. (2011). Geothermal Area Detection Using Landsat ETM+ Thermal Infrared Data and its Mechanistic Analysis- A Case Study in Tengchong, China. International Journal of Applied Earth Observation and Geoinformation, 13, 552-559. Ren, H., Du, C., Liu, R., Qin, Q., Guangjian, Y., & Li, Z. L. (2015). Atmospheric water vapor retrieval from Landsat 8 thermal infrared images. American Geophysical Union - Journal of Geophysical Research: Atmospheres, 120(5), 1723-1738. Rodríguez Rodríguez, G. F., & Vallejo Rodríguez, E. F. (2013). Sondeos Térmicos Superficiales en el Área Geotérmica de Paipa (Boyacá). Colombian Geological Survey. Research and Exploration of Geothermal Resources, Basic Geology, Bogotá, Colombia. Rodríguez, C. (1998). Exploración y Caracterización de los Recursos de Aguas Termominerales en el Municipio de Paipa. Corpoboyacá, Bogotá, Colombia. Rueda Gutíerrez, J. (2016). Cartografía de los Cuerpos Dómicos del Área Geotérmica de Paipa. Colombian Geological Survey. Research and Exploration of Geothermal Resources, Basic Geosciences, Bogotá, Colombia. Sekertekin, A., & Arslan, N. (2019). Monitoring Thermal Anomaly and Radiative Heat Flux Using Thermal Infrared Satellite Imagery – A case study at Tuzla Geothermal Region. Geothermics, 78, 243-254. Sobrino, J. A., Jiménez Muñoz, J. C., & Paolini, L. (2004). Land surface temperature retrieval from LANDSAT TM 5. Remote Sensing of Environment, 90, 434-440. United States Geological Survey (USGS) (2016). Landsat 8 (L8) Data Users Handbook. Technical Manual, South Dakota, United States. United States Geological Survey (USGS) (2018). Landsat 7 (L7) Data Users Handbook. Technical Manual, South Dakota, United States. Ussher, G., Harvey, C., Johnstone, R., & Anderson, E. (2000). Understanding the Resistivities Observed in Geothermal Systems, Proceedings World Geothermal Congress 2000, Kyushu-Tohoku, Japan, May- June, 1915-1920. Vásquez, L. (2002). Exploración Geofísica con el Método Electromagnético en el Municipio de Paipa. Colombian Geological Survey (Previously INGEOMINAS). Geo-Resources Branch, Geothermal Exploration Project. Bogotá, Colombia. Vásquez, L. (2012). Aplicación de Métodos Potenciales en el Área Geotérmica Paipa-Iza. Colombian Geological Survey (Previously INGEOMINAS). Subsoil Resources Branch. Bogotá, Colombia. Velandia, F. (2003). Cartografía Geológica y Estructural Sector Sur del Municipio de Paipa. Colombian Geological Survey (Previously INGEOMINAS). Geodynamic Project. Bogotá, Colombia. Wu, W., Zou, L., Shen, X., Lu, S., Su, N., Kong, F., & Dong, Y. (2012). Thermal infrared remote-sensing detection of thermal information associated with faults: A case study in Western Sichuan Basin, China. Journal of Asian Earth Sciences, 43(1), 110-117. Yu, X., Guo, X., & Wu, Z. (2014). Land Surface Temperature Retrieval from Landsat 8 TIRS—Comparison between Radiative Transfer Equation-Based Method, Split Window Algorithm, and Single Channel Method. Remote Sensing, 6(10), 9829-9852. Zhou, Y. R. (1998). The Application of Thermal Infrared Remote Sensing Techniques in Geothermal Surveying. Remote Sensing for Land & Resources, 67, 24-28.