Detection of geothermal potential based on land surface temperature derived from remotely sensed and in-situ data

ABSTRACT

Geothermal energy is a renewable and environmentally sustainable resource of increasing importance. However, areas with geothermal potential are not easily detected by traditional field investigations, requiring the development of new, robust, and reliable models for detection. In this study, remote sensing data and ground-based variables were used to detect and analyze geothermal resource potential areas. General Land Surface Temperature (GLST) was integrated using 5 years of remote sensing data. Landsat 8 daytime GLST (Landsat-GLST), Moderate Resolution Imaging Spectroradiometer (MODIS) daytime GLST (MODIS-DLST), and MODIS nighttime GLST (MODIS-NLST) data were integrated with Landsat Nighttime Land Surface Temperature (Night-LST), which not only filled the gap of Landsat Night-LST but also improved the spatial resolution of MODIS nighttime temperatures. Specifically, three independent variables (Night-LST, Distance From Known Geothermal Resource Points [DFGP], and Distance From Geological Faults [DFF]) were used to develop a weighted model to form a Geothermal Detection Index (GDI) based on Principal Component Analysis (PCA). Along with field verification, the GDI was successfully used to identify three geothermal activity areas in Tengchong City, Yunnan Province. Overall, this work provides a novel method for detecting geothermal potential to support the successful exploitation of geothermal resources.

KEYWORDS:

• Geothermal Detection Index (GDI)
• geothermal potential
• Principal Component Analysis (PCA)
• Land Surface Temperature (LST)

1. Introduction

With the recent increase in serious environmental problems, such as climate change and the depletion of traditional resources, the development of renewable energy sources, such as geothermal, solar, wind, and biomass, has intensified considerably. In particular, geothermal resources offer pollution-free and stable energy sources (T. Abbasi and S. Abbasi 2012; Kurek et al. 2021; Rezaie and Aghajani 2013; Van der Meer et al. 2014) that can effectively meet energy demands with minimal environmental impact (Chen and Xu 2013). Geothermal energy has been successfully harnessed on a large scale for space heating, industry, and electricity generation for many years (Ingvar 2001); however, uncertainty encountered in previous geothermal investigations has raised many problems (Jennejohn 2009). Therefore, it is particularly important to develop robust models that can detect geothermal activity and resource potential.

Previous studies have shown that suitable assessment criteria for geothermal potential include distance to mineral alteration zones, volcanic rock, faults, and hot springs; earthquake distributions; and magnetic anomaly distributions (Sun et al. 2012; Meng et al. 2021; Moghaddam et al. 2013, 2014; Mcguire et al. 2015; Li and Zhang 2017; Zhang et al. 2020). However, active geothermal fields often exhibit temperature anomalies at the ground surface because of heat and mass transfer processing from volcanic activity, hot springs, mud pots, or fumaroles (Tian et al. 2015). Currently, most studies detect geothermal activity by focusing on two criteria, namely geological structures and Land Surface Temperature (LST). Indeed, LST is a key parameter in the physics of land surface processes at regional and global scales (Zhang and He 2013; Wang et al. 2019; Gemitzi, Dalampakis, and Falalakis 2021; Li et al. 2021; Li, Xu, and Yao 2021). For example, Saepuloh et al. (2020) using The Thermal Infrared Red (TIR) of the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) data, Synthetic Aperture Radar (SAR) images, and ground data explored the relationship between surface temperature and roughness and identified the geothermal steam spots. Romaguera et al. (2018) explored the possibility of using two types of LST time series to detect geothermal anomalies and extract the geothermal component of LST.

The geothermal potential of areas within the Akarcay Basin in Turkey was also determined using a GIS-based Multi-Criteria Decision Analysis Method (MCDA), with hot springs reliably predicted as having geothermal potential (Yalcin and Gul 2017). Similarly, in Hokkaido, Japan, a geostatistical approach combined with TIR remote sensing data and three-dimensional (3D) temperature modeling based on well measurements was used for geothermal exploration (Tian et al. 2015). Zhang et al. (2020) identified and classified potential geothermal areas based on geothermal, geological, and geophysical spatial associations in Fujian Province, China. Based on the GIS supported weight-of-evidence model and fuzzy logic model, the geothermal prediction maps were established using the impact factor maps.

Rich data on the causes of Tengchong magma heat source geothermal resources. Qin et al. (2011) inversed the surface temperature by using the single-channel algorithm on the thermal infrared band, and identified four geothermal regions consistent with the fault development in Tengchong. Zhang et al. (2012) used remote sensing and Geographic Information System (GIS) for geothermal exploration in Tengchong, Yunnan Province, SW China. Xiong, Chen, and Huang (2016) used elevation information partition to identify the abnormal surface temperature areas and circled the geothermal abnormalities in Tengchong; Tang, Li, and Liu (2017) determined geothermal abnormal areas in Tengchong by eliminating the shadow influence and inversing the temperature.

In particular, Landsat datasets provide a unique opportunity for global-scale surface temperature assessment with more than 30 years of archived imagery, including TIR (Roy et al. 2014; Weng and Fu 2014; Weng, Fu, and Gao 2014; Cristóbal et al. 2018; Saepuloh et al. 2020). Over the past decade, LST derived from Landsat remote sensing data has been successfully used to assess geothermal potential (e.g. Li et al. 2014; Gemitzi, Dalampakis, and Falalakis 2021). More recently, LST methods including the split-window method (Becker and Li 1990; Wan and Dozier 1996), the single-channel method (Jiménez-Muñoz and Sobrino 2003), and radiative transfer equation-based methods have been developed. Moreover, obtaining LSTs using the radiative transfer equation-based method based on Band 10 has been verified as the approach providing the highest level of accuracy with Root Mean Square Errors (RMSEs) of <1 K (Yu, Guo, and Wu 2014). LST products of MODIS have been widely used in a range of research fields. This includes Li et al. (2012), who used multi-temporal nighttime MODIS LSTs to identify geothermal anomalies that show good agreement with relative geothermal gradient measurements. Furthermore, using historical data, LST models can be built to predict daytime temperature (Hengl et al. 2012). MODIS LST data have also been successfully combined with Landsat 8 data to improve the accuracy of LST measurements derived from satellites (e.g. Zhao et al. 2020).

In this study, we sought to build on this existing work by improving data quality and accuracy to better inform geothermal exploration. Specifically, as Landsat 8 LST data have a high spatial resolution, and MOD11A2 LST products contain both daytime and nighttime data, we combined these products to retrieve high-resolution nighttime LSTs. In addition, we developed a weighted principal component model to evaluate geothermal potential by combining Night-LST with geotectonic factors, which concluded that DFGP and DFF are useful indicators of geothermal potential. Tengchong City, Yunnan Province, located in Southwest China, is extremely rich in geothermal resources due to its geological structure. We take this county as an example for research.

2. Study area

Tengchong County is located in southwest Yunnan Province (Figure 1), Southwest (SW) China, with a total area of 5,845 km². This area sits at the intersection of the Eurasian and Indian Plates with 84% of the land area covered by mountains (Wei et al. 2003). The terrain is high in the north and low in the south, and the mountains run Northeast-west and North-south, which is basically the same with the regional structural line. The formation development in the study area was incomplete due to the new and Mesozoic magmatic activity, the lower Paleozoic formation was absent in the study area, only development Proterozoic rocks, Paleozoic rocks, Mesozoic rocks, Neogene rocks, and Quaternary sediments. Formation exposure is relatively sporadic and has poor continuity due to granite intrusion in different periods. Quaternary volcanic and fault structures are the most prominent characteristics of the region at present. The active faults in the Northeast, North, Northwest, and South of the region have developed well and are conducive to geothermal formation (Du et al. 2005）. The county is characterized by frequent tectonic-thermal events, seismic activity, and active faults, with more than 70 volcanic craters of different sizes and 140 hot springs (Wan, Zhao, and Kang 2005). As a famous volcanic and hot spring geothermal area, Tengchong displays a wide distribution of volcanic rocks and hot geothermal fields.

Figure 1. Location of Tengchong County, Yunnan Province, China.

3. Input ata

3.1. Remote sensing data

Remote sensing data, including Landsat 8 and MODIS, were the main focus of this study. The obtained images (Table 1) were free from cloud and snow cover. The geological structure map of the Tengchong area was obtained from the China Geological Survey, and land survey results for Tengchong County (2011) were obtained from the Land Management Department of Tengchong City. The verification data we used were from other studies (Li and Zhang 2017; Guo 2012; Ba 2017). Information about the Rehai geothermal field was also obtained from field surveys.

Table 1. Remote sensing data used in this study.

RS Data	Year	Month	Bands/Products	Spatial Resolution
Landsat 8	2015–2019	January/February	Band 4 (Red: 0.630–0.680 μm)/ Band 5 (NIR: 0.845–0.885 μm)/ Band 10 (TIRS 1: 10.60–11.19 μm)	30 m/100 m
MODIS	2015–2019	January/February	MOD11A2 Daily LST data	1 km

3.2. Distance from sites with known geothermal activities (DFGP)

To date, 64 sites of geothermal activity have been identified in the Tengchong area. These sites are mainly hot springs that have been explored by residents and opened to the public as tourism sites. Here, we determined 60 sites with strong geothermal activity, with a high probability that currently unknown geothermal activity exists around the known geothermal resource points. Therefore, the probability of geothermal activity decreases with distance from existing geothermal hotspots. Based on this, we generated a DFGP raster file showing the distance between each pixel and known geothermal hotspots (Figure 2(a)).

Figure 2. Raster files showing (a) distance from known geothermal resource points (DFGP) and (b) distance from faults (DFF).

3.3. Distance from faults (DFF)

Previous studies of geothermal resources in the Tengchong area have found that the fault structure controls the spatial distribution of geothermal activity, and high-temperature geothermal areas are concentrated at the intersection of faults. In the studied region, the active faults are predominantly north–south oriented (Shangguan, Bai, and Sun 2000). Correspondingly, the spatial distribution of geothermal heat is banded and curved following an anastomosing pattern consistent with the fault-line geomorphology. We created a DFF raster file based on the distance of each pixel from geological fault lines (Figure 2(b)). Note that fault depth was not considered.

4. Methodology

4.1. Technical framework

Night-LSTs were obtained by integrating Landsat 8 daytime General Land Surface Temperature (GLST; Landsat-GLST), MODIS Daytime GLST (MODIS-DLST), and MODIS Nighttime GLST (MODIS-NLST) with DFF and DFGP input data. From this, we developed a geothermal detection model based on a PCA weighting method (Figure 3).

Figure 3. Overall workflow adopted to detect geothermal potential.

4.2. Remote sensing data processing

4.2.1. LST retrieval

LST results from solar radiation and the energy emitted back from beneath the Earth’s surface. To detect geothermal resources, the effect of solar radiation can be subtracted from the total thermal radiation observed by the satellite sensor. The obtained thermal radiation intensity can then be converted into LST, which provides a scientific basis for the detection of potential geothermal resources. Current algorithms for retrieving surface temperature from remote sensing data include the radiative transfer equation method, single-window, and split-window algorithms, as well as multi-channel and multi-angle algorithms, among others. The radiative transfer equation method, also known as the atmospheric correction method, is mainly based on the composition of the thermal radiation intensity received by the thermal infrared satellite sensor. The algorithm subtracts the atmospheric influence factor and calculates surface temperature, using the relationship function between thermal radiation intensity and surface temperature. The radiation flux received by Landsat 8 Thermal infrared Sensor (TIRS) is mainly composed of surface heat radiation, atmospheric upward heat radiation, and atmospheric downward heat radiation, which are partially reflected by the sensor. The expression of the thermal infrared radiance received by the sensor can be written as (Sobrino, Jimenez-Munoz, and Paolini 2004):

(1) $L_{\lambda }=\tau \epsilon B\left(T_s\right)+\tau \left(1-\epsilon \right)L_{\downarrow }+L^{\uparrow }$(1)

where $L_{\lambda }$ is the radiance received by the thermal infrared sensor ($\mathrm{W}\cdot {\mathrm{m}}^{-2}\cdot \mathrm{s}{\mathrm{r}}^{-1}\cdot \mu {\mathrm{m}}^{-1}$), $\underset{\hspace{0.17em}}{\epsilon }$ is the specific surface emissivity, $B\left(T_s\right)$ is the thermal emission of a blackbody as expressed by Planck’s law, $L_{\downarrow }$ and $L^{\uparrow }$ are the downward and upward radiation intensities of the atmosphere, respectively ($\mathrm{W}\cdot {\mathrm{m}}^{-2}\cdot \mathrm{s}{\mathrm{r}}^{-1}\cdot \mu {\mathrm{m}}^{-1}$), and τ is the atmospheric infrared transmittance. The upward and downward radiation and transmittance of the atmosphere were obtained, using the atmospheric correction parameter calculator provided by the National Aeronautics and Space Administration (NASA) https://atmcorr.gsfc.nasa.gov (Barsi et al. 2005). Black body luminance can be obtained from the formula (2):

(2) $B{\left(T_s\right)}^1=\frac{\left[L_{\lambda }-\tau \left(1-\epsilon \right)L_{\downarrow }-L^{\uparrow }\right]}{\tau \epsilon }$(2)

The relationship between radiation energy, temperature, and wavelength within the specifications of TIRS 10 can be simplified according to Planck’s law, as follows:

(3) $B\left(T_s\right)=K_1/(e^{K_2/T_s}-1)$(3)

(4) $T_s=K_2/ln\left(K_1/B\left(T_s\right)+1\right)$(4)

where $T_s$ is the surface temperature, ${\mathrm{K}}_1=774.89\mathrm{W}\cdot {\mathrm{m}}^{-2}\cdot \mathrm{s}{\mathrm{r}}^{-1}\cdot \mu {\mathrm{m}}^{-1}$, $K_2=1321.08K$, based on previous studies (Van der Meer et al. 2014) on the range and spectral response function of TIRS 10.

Subsequently, specific surface emissivity $\epsilon$ and the radiance received by thermal infrared sensor $L_{\lambda }$ can be calculated. $L_{\lambda }$ can be directly obtained by radiometric calibration of the thermal infrared band 10. The processes of temperature retrieval by the radiation transmission equation method are shown in Figure 4. Land-specific surface emissivity, or emissivity, refers to the ratio of the radiation emitted by a land surface to the radiation emitted by a blackbody at the same temperature (Li et al. 2013).

Figure 4. The workflow of LST retrieval, using the radiative transfer equation-based method.

This ratio is determined by surface composition, roughness, the observed range of wavelengths, and other factors. Remote sensing images are mainly composed of water, urban land use, and natural land surfaces. Natural land surfaces, covered by different types of vegetation, differ in their specific emissivity. One of the most common classification indices, which was also used in the present study, is the Normalized Difference Vegetation Index (NDVI). As NDVI is a normalized index, the calculation error from atmospheric influence on NDVI is not very large. In the case of the present study, the surface temperature error caused by NDVI remained less than 0.1 K. NDVI can be obtained directly by the Operational Land Imager (OLI) B4 and B5 after radiometric calibration. The proportion of ground covered by vegetation is closely related to the temperature ratio and specific surface emissivity of land surfaces. Thus, NDVI can be used to estimate vegetation coverage $P_v$, as follows (Sobrino, Jimenez-Munoz, and Paolini 2004):

(5) $P_v=\frac{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}$(5)

where $\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ represents the minimum value of NDVI without vegetation, and $\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum value for an area completely covered by vegetation. For the calculated NDVI values, a certain amount of noise is inevitable, for this reason $\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ and $\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ generally are reported within their respective confidence ranges. However, for the calculation of vegetation cover using equation (5)(5) $P_v=\frac{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}$(5) , it is assumed that the pixel value is 1 if NDVI>$\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ and it is 0 if NDVI<$\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$.

For the calculation of land surface emissivity, the specific emissivity of typical ground objects can be represented by an estimate, and the band 10 thermal infrared band of TIRS and Thematic Mapper/Enhanced Thematic Mapper (TM/ETM +) 6 thermal infrared band have similar spectral range, so the same calculation method of surface emissivity of TM/ETM+6 is used in this study. In addition, Yu, Guo, and Wu (2014) used a single band for surface temperature inversion, and found that the band 10 was relatively Band 11 has high inversion accuracy.

In this study, LST products of MODIS were used to evaluate if the LST results from Landsat imagery were reliable and used to synthesize the night land surface temperature data combining the Landsat land surface temperature. MOD11A2 Daytime LST products provide per-pixel Land Surface Temperature and Emissivity with 1-km spatial resolution. It included Daytime Land Surface Temperature, Nighttime Land Surface Temperature, Daytime LST Quality Indicators, and other data layers. To obtain LST data, it is necessary to extract the corresponding layers from MOD11A2 products. For this study, the MODIS Reprojection Tool (MRT) (Dwyer and Schmidt 2006) software was used to re-project the MOD11A2 data and convert the default Hierarchy Data Format-Earth Observation System (HDF-EOS) data format to GeoTIFF. To reduce the required storage space, the LST product needed to undergo a numerical conversion process. To obtain LST in degree Celsius, MOD11A2 raw data had to be converted as follows:

(6) $T_m=DN\times Scalefactor-273.15$(6)

where $T_m$ is the ground temperature in degree Celsius, DN is the raw value from the LST_Day_1km layer of the MOD11A2 product, and $Scalefactor$ is the proportional coefficient for the product, in this case, 0.02.

MODIS LST products were used to evaluate whether the LST results obtained from the Landsat imagery were reliable, which were then used to synthesize the Night-LST data. First, the corresponding layer was extracted from the multiple data layers of the MOD11A2 LST product. The MRT was then applied to convert the default data format to GeoTIFF. The MOD11A2 raw data were also converted to degree Celsius.

4.2.2. Construction of GLST

As the LST data for 1 year cannot accurately reflect true LST due to complex variability over a range of timescales, a composite of remote sensing data for 5 years was used to synthesis based on the standard deviation. Thus, the standard deviations for LST from 2015 to 2019 were used to determine surface temperature anomalies at different points in time, thereby representing more generalized surface temperature conditions. The weighted and synthesized LST estimates were calculated as follows:

(7) $W_i=\frac{S_i}{S_1+S_2+\dots +S_i}$(7)

(8) $GLST=W_1{T}_1+W_2{T}_2+\dots +W_i{T}_i$(8)

where $W_i$ is the weight of the ith image, $S_i$ is the standard deviation of the ith image, $T_i$ is the LST obtained from the i-th image, and GLST is the comprehensive temperature.

Comparisons of Landsat-GLST, MODIS-DLST, and MODIS-NLST data are shown in Figure 5 and Table 2. Although the selected data were from the same period, there was some deviation due to slight differences in the acquisition times and temperature-retrieval algorithms. Nevertheless, in general, the spatial pattern of pixel values was in good agreement between Landsat-GLST and MODIS-DLST. High-temperature anomalies were concentrated in similar regions, indicating that temperature results are reliable. Therefore, the spatial distribution of the measured temperature range still has reference value, and the comparison demonstrated that the results from remote sensing data were reliable, especially the occurrence of high-temperature anomalies in the eastern farmland area. For MODIS-NLST, the spatial pattern was different from that of the DLST values, which reflects the intensity of insolation, the different types of land cover, and the influence of human activities. Take Rehai geothermal field as an example in Figure 5, the Landsat-GLSTs and MODIS-DLSTs are broadly similar, although there are still some differences due to the higher spatial resolution of the MODIS-DLST dataset. In comparison, the pattern of MODIS-NLSTs is quite different from the other two datasets. This shows that the surface temperature during the day is affected by other external factors and cannot accurately reflect the geothermal conditions, so it once again proves that the surface temperature at night is more reflective of the abnormal geothermal activities. Compared to the MODIS-DLSTs in Table 2, the temperature range is greater in the Landsat dataset, with a maximum of 30.34°C and a minimum of 5.29°C, although the average temperature is similar to that of the MODIS-DLST dataset. As expected, the nighttime LSTs are generally much lower than the daytime values. As a result, these statistical data once again prove that the use of Landsat data can more effectively detect geothermal abnormal areas.

Figure 5. Comparison of GLST (a) Landsat-GLST; (b) MODIS-DLST; and (c) MODIS-NLST.

Table 2. Summary statistics for Landsat-GLST, MODIS-DLST, and MODIS-NLST datasets.

GLST Data	Maximum (℃)	Minimum (℃)	Mean (℃)		Standard Deviation
Landsat-GLST	30.34	−5.29		13.24	3.39
MODIS-DLST	23.45	4.67		14.64	2.83
MODIS-NLST	9.62	−2.19		4.46	1.89

4.2.3. GLST verification

To verify the LST estimates, the RMSE (Fridleifsson 2001), Mean Absolute Error (MAE) (Willmott and Matsuura 2005), and Structural Similarity Index Measure (SSIM) (Qin et al. 2004) were used to evaluate the consistency between the GLST values. For the comparison of Landsat-GLST and MODIS-DLST, the RMSE indicated abnormal values in the image. In comparison, the MAE was less than 3°, indicating that the errors of the two images were small. Finally, the SSIM was 0.91, which shows that the structural similarity between the two images was very high (Figure 6(a)). Based on these results, the reliability of the temperature retrieval method was confirmed, including support its use for the synthesis of Night-LSTs. In addition, the difference between MODIS daytime temperature data and nighttime temperature data was more pronounced in high-temperature regions.

Figure 6. (A) Comparison of Landsat-GLST and MODIS-DLST; (b) Comparison of MODIS-DLST and MODIS-NLST.

4.2.4. Construction of Landsat Night-LST

Based on the comparison and verification of the GLST data and to eliminate the influence of the spatial resolution of the MODIS dataset and compensate for the lack of Landsat nighttime surface temperature data, Landsat 8 nighttime LST inversion was performed as follows:

(9) $W_m=\frac{M_n-M_d}{M_d}$(9)

(10) $L_n=L_d+W_m\times L_d$(10)

where $W_m$ is the weight, $M_n$ is MODIS-NLST, $M_d\hspace{0.17em}$ is MODIS-DLST, $L_n$ is the Night-LST, and ${\mathrm{L}}_{\mathrm{d}}$ is Landsat-GLST.

4.3. PCA

PCA is a powerful multivariate statistical method that is commonly applied in a wide range of areas (Bro and Smilde 2014). PCA can reassemble a group of multiple independent variables into a new, smaller set of variables that retain as much of the information contained in the original data as possible. The new, condensed variables are called principal components. In the process of synthesizing multiple independent variables to detect geothermal resources, it is difficult to determine the magnitude of each factor’s role in the process, which affects the detection capability of final results. PCA can effectively eliminate redundant information and determine the impact weight of each influencing factor to construct the detection model. Before performing PCA, an original sample matrix needs to be constructed. For this purpose, the study area was randomly divided into n blocks. A matrix X is established by using the average values of N influence factors in n regions as an index for the principal component analysis. The ith column vector of matrix X is the sample column vector $x_i$ (i = 1,2, … N). The column vectors are used to compute the correlation coefficient matrix R of matrix X as follows:

(11) $R_{N\times N}=\frac{Cov\left(X_i,X_j\right)}{\sqrt{D(X_i)}\sqrt{D(X_j)}}\hspace{0.17em}(i,j=1,2\dots ,N)\hspace{0.17em}$(11)

where $\mathrm{C}\mathrm{o}\mathrm{v}\left(X_i{X}_j\right)$ is the covariance matrix, composed of the covariances of $X_i$ and $X_j$. $D\left(X_i\right)$ and $D\left(X_j\right)$ are the variances of the corresponding column vectors. The column vectors can also be normalized to mean 0 and variance 1, resulting in a correlation coefficient matrix:

(12) $X_i^{\prime }=\frac{X_i-Mean\left(X_i\right)}{Std\left(X_i\right)}$(12)

(13) $R_{N\times N}=Cov\left(X_i^{\prime },X_j^{\prime }\right)$(13)

where $\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}\left(X_i\right)$ and $Std\left(X_i\right)$ are the averages and standard deviations of the sample variables $X_i$, and $X_i^{\prime }$ is defined as a positive standard matrix. Subsequently, the eigenvalues and eigenvectors are obtained from the correlation coefficient matrix R by solving the equation:

(14) $\left|R-\lambda E\right|=0$(14)

(15) $\left(R-{\lambda }_{i}E\right)x=0$(15)

where $\lambda$ is the N eigenvalues of the matrix R and x is the eigenvector that corresponds to the eigenvalue ${\lambda }_i$(i = 1,2, …, N). Arranging the eigenvalues from largest to smallest (${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_N$), the variance contribution rate g and the cumulative variance contribution rate G are calculated. The higher the cumulative variance contribution rate is, the lower the number of resulting principal components, which in turn indicates that less information has been lost. In general, eigenvalues with a cumulative variance contribution over 80% should be selected as main components, as in this case the quantity m of the principal component should be satisfied:

(16) $\frac{{\sum }_{i=1}^m{\lambda }_i}{{\sum }_{i=1}^N{\lambda }_i}=\frac{g}{G}\ge 85\mathrm{\%}\left(m\le N\right)$(16)

The eigenvector $x_i$ corresponding to m principal components is used as the coefficient of the comprehensive score model, and the weight of each influencing factor can be determined as follows:

(17) $W_i=\frac{{\sum }_{j=1}^m{x}_{ij}{g}_j}{{\sum }_{j=1}^m{g}_j}$(17)

(18) $w_i=\frac{W_i}{{\sum }_{j=1}^N{W}_j}$(18)

where $x_{ij}$ is the linear combination coefficient of the i-th influencing factor in the j-th principal component, $g_j\hspace{0.17em}$ is the variance contribution rate of the j-th principal component, $w_i$ is the standardized result of the weight $\hspace{0.17em}{W}_i$ of the i-th influencing factor.

4.4. GDI

The GDI was calculated based on the weight of each influencing factor, where the weighted average operator of each factor is:

(19) $GPI={\sum }_{i=1}^N{w}_i{c}_i$(19)

where $c_i$ is an influencing factor for detecting geothermal activity, and $\hspace{0.17em}{w}_i$ represents the weight of each factor. Based on the PCA results, DFF, DFGP, and Night-LST were selected as the influencing factors to calculate the GDI.

5. Results and validation

5.1. Spatial pattern of Night-LST

Night-LST data are shown in Figure 7(a) and summarized in Table 3. Compared with the MODIS-NLSTs, Night-LST shows good spatial consistency with similar high-temperature anomalies in association with some human activity and farmland areas in the south. The temperature range of the Night-LST dataset is relatively large, but the average and standard deviation are broadly the same as the MODIS-NLST dataset. There is a certain degree of difference in the spatial pattern of LSTs between the nighttime and daytime values. This difference may be due to the impact of human activity and the heat island effect, i.e. the relatively higher temperature of urban surfaces in some cities compared to other areas. Thus, the inversion data of the Night-LST were available. Importantly, based on the identification of areas with abnormal Night-LSTs, geothermal resource potential can be effectively detected.

Figure 7. (A) Spatial pattern of Night-LST; (b) Classification of land cover based on Random Forest classification. The two images in the middle are detailed enlargements of the areas indicated.

Table 3. Summary statistics for Night-LST.

Data	Maximum (℃)	Minimum (℃)	Mean (℃)	Standard Deviation (℃)
Night-LST	14.25	−6.21	4.04	1.83

Figure 7(b) shows a land-use classification map based on Landsat 8 multispectral data combining the NDVI, Normalized Difference Built-up Index (NDBI), texture analysis, and regions of interest obtained using the Random Forest classification. Figure 7(b) is used as a qualitative analysis to assist in interpreting the types of land cover in the study area. Based on the Night-LST dataset, the high-temperature areas are mainly concentrated in the vegetated and farmland areas as well as around some cities in the south. For example, the Panzhihua geothermal field in the south with high-temperature anomalies (Figure 7) is mostly by natural vegetation and cropland.

5.2. The results of PCA

I propose creating a map of 10 × 10 fishing net statistical zones, generating 78 regions, and a correlation coefficient matrix was established for the independent variables. The correlation coefficients and eigenvalues along with their variance contribute rates are shown in Table 4, and eigenvectors are shown in Table 5. The resulting two eigenvalues had a cumulative variance contribution rate of 88.6%. The weights of the first two component influencing factors were calculated using their corresponding two eigenvectors (Table 5). According to previous studies, DFF, DFGP, and Night-LST have associations with geothermal resources, and they are similar to the content expressed by the first principal component. Therefore, the first principal component was considered to contain information relating to three influencing factors.

Table 4. Eigenvalues of principal components and their variance contribution rates.

Component	Eigenvalue	Percentage of Variance (%)	Cumulative Percent (%)
1	1.737	57.9	57.9
2	0.920	30.7	88.6
3	0.343	11.4	100.0

Table 5. Weights of the influencing factors and coefficients of the principal components.

Evaluation Factor	Coefficient of Component 1	Coefficient of Component 2
DFF	0.32	0.95
DFGP	0.66	−0.27
Night-LST	0.68	−0.18

5.3. Spatial pattern of GDI values

The GDI results are shown in Figure 8, with high values concentrated around cropland and urban areas, which indicates that (i) human activity is concentrated in areas of higher geothermal activity and (ii) the heat produced by human activity may alter Night-LST patterns. Outside of the urban areas, high GDI values indicate key areas for further geothermal exploration, with three geothermal potential areas identified. The first area (area “a” in Figure 8) is centered on the Langpu and Rehai thermal fields, mainly covered by vegetation and with human activity, with existing developed geothermal activity points. The second area is located around the Panzhihua and Heishihe thermal fields (area “b” in Figure 8), where geothermal activities are already developed and faults are relatively concentrated. The third area is located around Beihai thermal field (area “c” in Figure 8) in association with the north-south Dazhai-Longjiaying-Jiuzhai Fault. There are also some areas with the potential for geothermal exploitation in the Ruidian thermal field.

Figure 8. Spatial distribution of Geothermal Detection Index (GDI) values. The enlargements show areas marked “a”, “b”, and “c” in the left image centered on thermal fields.

Based on Figure 8, the GDI values in Tengchong range between 0.3 and 1.6 and the area with values between 1.3 and 1.6 were chosen as the area of geothermal development potential. Land-use data for this area are shown in Table 6, with almost half of the pixels indicating vegetation cover and approximately 41% cropland cover; areas of other human activity accounted for less than 10%. The altitude range of Tengchong County is 971.76 m to 3671.19 m. Compared to the elevation of the entire county, the average elevations of areas with GDI values ranging from 1.3 to 1.6 are relatively close, with the highest average elevation located in vegetation areas, followed by cropland. The lowest average elevation is in waterbody, followed by buildings. This is more in line with reality. For the average GDI value, the five types of land cover are 1.36 and 1.37. Therefore, the data analysis in the comprehensive table, most of these areas, which indicate rich geothermal resources, has not yet been developed. In view of some vegetation-covered areas with high altitude are difficult to develop, and the large-scale development in areas with frequent human activities is temporarily not carried out, the remaining vegetation-covered areas, cropland, bare soil, and some waterbodies have great potential for geothermal development.

Table 6. Land-use classifications and corresponding pixel counts, proportion (%), area (m²), average altitude (m) and average GDI within GDI values ranging 1.3–1.6.

Land-Use Classification	Pixel Count	Percentage (%)	Area (m^2)	Average Altitude (m)	Average GDI
Vegetation	113445	49.8	102100500	1768.38	1.36
Cropland	94418	41.4	84976200	1624.20	1.37
Waterbody	452	0.2	406800	1234.11	1.36
Buildings	7483	3.3	6734700	1441.99	1.36
Nudation	12132	5.3	10918800	1553.90	1.37

5.4. Validation

5.4.1. Langpu geothermal field

The distribution of GDI values in the Langpu geothermal area is shown in Figure 9(a), ranging from 1.2 to 1.4. Using these values, we tested three thresholds (1.34, 1.36, and 1.38) for the extraction of geothermal areas with the greatest development potential (Figure 9). The area corresponding to a threshold of 1.36 (Figure 9(c)) was most suitable, producing an area that would neither require excessive survey work (as is the case for the 1.34 threshold, Figure 9(b)) nor exclude potentially valuable locations (as is the case for the 1.38 threshold, Figure 9(d)). This area contained 654 favorable points with a total area of 0.59 km². Compared to the recorded first-grade geothermal area of 0.6 km² (Zhang et al. 2012), the extraction accuracy was 98%.

Figure 9. GDI values in the Langpu geothermal field (a) and the geothermal potential detection area extracted using three thresholds of 1.34 (b), 1.36 (c), and 1.38 (d).

5.4.2. Rehai geothermal field

Temperature was measured in eight drill wells in the Rehai area (Figure 10), and the corresponding GDI values are shown in Table 6. In general, the GDI values (ranging from 1.317 to 1.378) increased with an increase in well temperature (ranging from 30°C to 96°C), but there were some exceptions, which may reflect differences in the well depths (Table 7). The locations of the eight measured drilling points are shown in Figure 11(a), with Hamazui having the highest temperature (up to 96°C) and GDI value (1.378). Figure 11(b) shows GDI values around Da Zhai with a temperature of 71°C. Figure 11(c) shows GDI values around Shangqiluo, with a temperature of 30°C, indicating that the GDI results were also reliable for detecting potential geothermal areas with relatively low temperatures.

Figure 10. Field surveying in the Rehab area.

Figure 11. (A) Locations of eight drilling sites and Geothermal Detection Index (GDI) values in the Rehai area; (b) GDI values around Dazhai (CK2); (c) GDI values around Shangqiluo (CK8).

Table 7. Average GDI values corresponding to drilling well temperatures.

ID	Drill Name	Well Temperature (℃)	GDI
CK1	Hamazui	96	1.378
CK2	Da Zhai	71	1.350
CK3	Shangjia Zhai 1 Well	61	1.376
CK4	Qiluoyingfang Well	56	1.364
CK5	Shangjia Zhai 2 Well	49	1.354
CK6	Langpu Zhai	45	1.350
CK7	Bangeshan	32	1.314
CK8	Shangqiluo	30	1.317

6. Discussion

6.1. Impact of human activity on LST retrieval

Interference from human activity can complicate the interpretation of LST results. Indeed, some important temperature anomalies indicating the presence of geothermal resources may be neglected if areas of human activity are simply excluded. For this reason, it is important to further consider land use in the study area. The 2011 Tengchong City land survey identified areas affected by intense human activity, as shown in Figure 12. A comparison with Figure 7(b) indicates that the survey results can be considered fairly accurate. In the future, these land-use data could be combined with visual interpretation to determine the land-cover type of study sites as training samples for machine learning, which could allow the more effective exclusion of human interference and target these sites for subsequent research.

Figure 12. Areas of intensive human activity in the Tengchong region.

In addition to this, the urban heat island effect also needs to be taken into account. Urban Heat Islands (UHI) are a manifestation of these changes, a phenomenon in which urban areas are warmer than surrounding rural areas (Mirzaei and Haghighat 2010). Urbanization creates a process in which natural landscapes are replaced by built-up areas. Therefore, it affects the absorption, storage, and emission rates of surface energy and releases excess anthropogenic energy (Hart and Sailor 2009). Urban heat islands are one of the main negative impacts of increased human activity and impervious surface cover in urban environments (Firozjaei et al. 2020; Huang, Liu, and Li 2021; Jiang et al. 2022). As an important part of the urban impervious surface, buildings and residences have different surface radiation and energy distribution from the natural surface. In order to ensure the accuracy of the identification of the geothermal potential area by the LST. In the future, it can be combined with Night-Time Light data or other effective methods to identify and exclude as much as possible, which is still a complex and challenging task (Liu, Chi, and Kuang 2014; Peng et al. 2016; Kuang 2018; Shao, Wu, and Li 2021; Qian et al. 2022; Yu et al. 2022; Li et al. 2022; ElGharbawi and Tamura 2022).

6.2. Influence of topography and fault depth on geothermal detection

LST is controlled by the complex interplay of topography, incoming radiation, and atmospheric processes as well as by soil moisture, land cover, and vegetation types, especially in mountainous areas (Bertoldi et al. 2010; Jamali et al. 2022; Hussain et al. 2022). Therefore, to improve the reliability of remotely sensed LST data, topographic correction is needed to adjust the spectral characteristics of land surfaces and account for the differences in surface irradiation caused by terrain. In this context, shading also had important effects on the land-surface energy balance in mountainous terrain (Bellasio et al. 2005). Therefore, the influences of geological factors complicate the quantification of topographic effects on LST and, thus, limit model detection accuracy. Furthermore, faults with different activity depths result in variable surface temperatures and degrees of geothermal fluid activity. Deep faults can induce upwelling of geothermal water, while shallow faults only provide channels for lower-temperature groundwater. The burial depth of geothermal resources and, hence, the economic feasibility of exploitation efforts. Therefore, further research is required to examine the influence of fault depth on geothermal potential detections.

6.3. Influence of land-surface thermal properties

Energy at the land surface is mainly derived from solar radiation, bioenergy, and the internal heat of the Earth (Zhang and Cao 2003). As we specifically aimed to explore the relationships between heat from the internal Earth and land-surface temperature, we selected MODIS-NLST products to reduce the influence of daytime heat sources. However, given the range of land-cover types both in our study area and elsewhere, the thermal properties of the land surface can vary greatly (Zhang et al. 2015). Therefore, differences in the thermal capacity, conductivity, and diffusivity of different land-cover types should be taken into consideration in future research.

7. Conclusions

We developed and applied a new method for detecting geothermal activity using Tengchong as a case study. We derived high-resolution Night-LST data with satisfactory verification results. These data were combined with DFF and DFGP datasets in a multi-factor weighted model based on PCA, yielding a comprehensive GDI for assessing geothermal resources. By analyzing the spatial distribution of GDI values, we identified three geothermal potential areas in Tengchong, including vegetated areas that have not yet been developed. These locations represent potentially valuable areas for the development and utilization of geothermal energy resources. Considering that areas with high LST coincided with areas of high human activity, additional research is now needed to examine human land uses on the derived GDI results. In addition, a greater range of variables linked to geothermal resources could be added to the model to improve accuracy and detection capability.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The remote sensing data are available at [http://glovis.usgs.gov/], the fault data and the sites of geothermal activity are derived from the following resources available in the public domain: [https://data.earthquake.cn/index.html, https://geocloud.cgs.gov.cn/]. The in-situ data that support the findings of this study are available from the corresponding author [Zhifang Zhao], upon reasonable request.

Additional information

Funding

This research is supported by the National Natural Science Foundation of China [Grant No. 41961064]; Yunnan Fundamental Research Projects [Grant No. 202001BB050030]; the Plateau Mountain Ecology and Earth’s Environment Discipline Construction Project [Grant No. C1762101030017]; Joint Foundation Project between Yunnan Science and Technology Department and Yunnan University [Grant No. C176240210019]. Joint Foundation Project between Yunnan Science and Technology Department and Yunnan University [Grant No. 2018FY-019].

Notes on contributors

Fei Zhao

Fei Zhao is an associate professor in the School of Earth Sciences, Yunnan University, Kunming, China. He received his PhD degree from Wuhan University, Wuhan, China, in 2012. His main research interests include cartography, spatial analysis, and urban sustainable development.

Zhiyan Peng

Zhiyan Peng is a postgraduate in School of Earth Sciences, Yunnan University, Kunming, China. Her research interests are remote sensing application and spatial analysis.

Jiangkang Qian

Jiangkang Qian is a PhD student in Cartography and GIS at Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing, China. His research interests are the quantification of anthropogenic heat and urban thermal environment simulation.

Chen Chu

Chen Chu is a master student majoring in Geographical Information Sciences at Institution of Geographic Sciences and Natural Resources Research, Chinese Academy of Science, Beijing, China. His research interest is spatio-temporal data mining.

Zhifang Zhao

Zhifang Zhao is a professor in the School of Earth Sciences, Yunnan University, Kunming, China. Her research interests include resource and environmental remote sensing monitoring, mineralized alteration remote sensing anomaly enhancement extraction, and remote sensing monitoring of land and resources in border areas.

Jiangqin Chao

Zhifang Zhao is a professor in the School of Earth Sciences, Yunnan University, Kunming, China. Her research interests include resource and environmental remote sensing monitoring, mineralized alteration remote sensing anomaly enhancement extraction, and remote sensing monitoring of land and resources in border areas.

Jiangqin Chao is a Lecturer in School of Geographical Sciences and Tourism, Zhaotong University, Zhaotong, China. She is pursuing the PhD in Cartography and GIS at the Institute of International Rivers and Eco-Security, Yunnan University, Kunming, China. Her research interests are thermal infrared quantitative remote sensing inversion and thermal infrared geothermal geological exploration.

Shiguang Xu

Shiguang Xu is a professor in Yunnan Geology and Mineral Explorating & Engineering Corp. (Group) and a chief engineer in Yunnan Geological Engineering Survey and Design Institute. His research interests are water resources, geological environment, and prevention and control of geological disasters.