Retrieval of Land Surface Temperature from Landsat 8 Data of the Dandong-Liaoyang Geothermal Area

: In recent years, with the intensification of problems like the depletion of traditional fossil fuels and environmental degradation, the development of new energy sources has become a key long-term strategy in China. Geothermal energy has attracted much attention due to its advantages of abundance and low environmental impact. Based on infrared data sensed remotely by the Landsat 8 satellite, this paper reports a verification of the atmospheric-correction method for extracting the surface temperature of the Dandong-Liaoyang geothermal region in all months of 2014. The method combines the abnormal points in the inversion results with the local sites of hot springs, structures, local historical air temperatures, and land surface temperature/emissivity data (MOD11_L2). Results showed that these data sources were spatially distributed in similar ways, which indicates that these results can be used to identify promising geothermal resources from publically available thermal-infrared remote-sensing data.


Introduction
In recent years, as fossil fuels are growing scarce, researchers have become focused on finding more-reliable and sustainable energy resources. Geothermal energy is both reliable and sustainable and does not emit pollutants. The economically rational exploitation of geothermal resources can address the social challenges of multiple resource constraints and heavy environmental pressure [1-2]. Land surface temperature (LST) is a very important indicator in the exploration for geothermal resources. By identifying differences in thermal remote sensing data, researchers can identify the location and quantities of available geothermal resources. To obtain the spatial and temporal distribution of the surface temperature over the whole earth range or large regions, only thermal-infrared remote sensing data can be used. Thermal-infrared remote sensing data is collected from satellites or aircraft with sensors that use bands limited to infrared wavelengths emitted from objects on the surface. This data is correlated with objects and can then be used to infer the surface 2 temperature. At present, many satellites in orbit carry thermal-infrared sensors, such as AVHRR, ASTER, MODIS, Landsat TM/ETM+, and Landsat 8 thermal-infrared sensor (TIRS).
The rapid development of thermal-infrared geothermal retrieval technology allows the quantitative study of thermal resources and monitoring of changes in the distribution of small thermal fields. These data sources are characterized by high speed, low cost, and wide range. This breadth of data and analysis techniques has laid a solid technical foundation for the application of remote sensing technology to the exploration of geothermal resources. In recent years, several methods to infer surface temperature using thermal-infrared remote sensing data have come to prominence: the split-window algorithm [3]; the mono-window algorithm [4-6]; the multiple-channels algorithm [7][8]; the atmospheric-correction method (thermal channel method) [9-10].
The Dandong-Liaoyang region is in the east of Liaoning province, China, which is rich in geothermal resources and has a well-developed infrastructure. More than 20 hot springs are known in this area, including the famous Dongtang and Wulongbei springs, so the area has good contrast for testing the analysis of thermal-infrared remote sensing data. The thermal-infrared remote sensing data is selected by Landsat 8 data. The Landsat 8 satellite has an advanced TIRS with a spatial resolution of 30 m. Below, we report how we used the Landsat 8 thermal-infrared remote sensing data to compare the anomalous areas in the extracted surface temperature map with the spatial distribution of fault structures and known hot spring locations in the study area. This comparison gives evidence for the effectiveness of the temperature-retrieval method. Below, the temperature retrieval results are also compared with LST data from the National Aeronautics and Space Administration's (NASA) MODIS LST and Emissivity dataset (MOD11-L2), the distribution of fault zones and hot springs.
The rest of this paper is organized as follows: the section 2 introduces the geographic, geological, and geothermal characteristics of the study area; section 3 details our method for extracting LST from Landsat 8 data; section 4 discusses our comparative method to verify the accuracy of the retrieval algorithm; and section 5 concludes the paper.

Geography
The study area is the eastern part of Liaoning Province, China ( The geomorphology of the survey area is complex and mostly mountainous. The north part belongs to the hilly region of eastern Liaoning, which is the western extension of the Changbai Mountain range. The vegetation is lush, especially in summer. The topography is complex and varied, and the height difference between peaks is also very large. Most of the peaks are between 500-700 meters above sea level, and the highest peak is about 1180 meters above sea level. The south of this region features low hills and a coastal plain, with an elevation below 500 meters. The main mineral resources are gold, silver, lead, zinc, boron, iron boride, magnesite, iron, copper, pyrite, talc, jade, and limestone in more than 18 forms.
The crust in the study area is a binary structure with distinct base and cap, and its crystalline basement is widely exposed in the area, which is mainly composed of the Archean metamorphic pluton and Paleoproterozoic Liaohe groups. Caprock is exposed sporadically, mainly including Neoproterozoic, Paleozoic terrestrial clastic rocks, carbonate rocks, Mesozoic river and lake facies, and continental volcano-sedimentary series [12]. The study area, which is still in the active stage after many tectonic evolutions, features three large structure faults trending in the NNE, NE, and NW directions.

Geothermal characteristics
Due to the development of branching faults and the breakage of the crust structure, surface water can easily infiltrate deep underground, where it is heated and then rises to the surface to form hot springs. As Figure 3 shows, the study area in the eastern region has many hot spots. These are known to have several large hot springs, with spring water temperatures between 30 ºC -71ºC, which are relatively medium or low temperatures for hot springs. They are generally distributed in the NNE and NE of the Liaodong region. Geothermal activity is driven by major faults; most geothermal fields are located at the intersection of fault structures. Many geothermal fields have been found in the study area, most of which are exposed in the form of hot springs. These hot springs are mainly fed by bedrock fissure water and are distributed mainly in the mountainous uplift area. According to previous data, the NNE   In sum, the study area includes several active tectonic belts with faults, which provide channels for the formation of hot springs, and a large amount of Mesozoic granite rock serves as thermal storage media that provide stable and reliable sources of heat. These previously identified resources suggest that the study area is ripe for the development of geothermal energy plants [13].

Methods
Since the study area has four distinct seasons, with severe vegetation cover in summer and severe snow cover in winter, a suitable time and phase should be chosen when considering the selection of remote-sensing data for our analysis. Landsat 8 thermal-infrared remote sensing data from each month of 2014 was selected in this paper for reasons of accuracy, acquisition speed, and clear boundaries. After a review of previously discussed surface-temperature retrieval algorithms, we chose the atmospheric-correction method (also known as the heat-radiation-transfer equation algorithm) to infer surface temperature.
ArcGIS 10.1 and ENVI 5.1 software were used for comprehensive data analysis. A basic flow chart of our analysis is shown in Figure 4: To use the spectral information from the remote sensing images effectively, the radiation intensities need to be calibrated. For the Landsat 8 data, the ENVI software includes an automatic calibration tool. Since the purpose of radiometric calibration in geothermal retrieval is to convert DN value to radiance values, we opened the radiometric calibration settings in the toolbox and changed the calibration type to radiance. When other default parameters are selected, the tool will deliver Band10 radiance images. If the Landsat 8 is marked with L1T when downloaded, it indicates that the data has undergone geometric correction with topographic features at the time of acquisition, and can generally be used directly.

LST retrieval
Using the basic principles of radiative transfer equations for the temperature retrieval, an atmospheric model is used to simulate the relevant influences of the atmosphere. First, the algorithm estimates the atmosphere's effects on thermal radiation or absorption transfer at the surface. Next, that part of the total thermal radiation detected by the sensor is eliminated, then the thermal radiation intensity from the surface is calculated, and finally, this intensity is converted to temperatures using the ENVI software kit [17][18]. The details of this process are given in a flowchart in Figure 5. As radiation is transmitted from the earth through atmospheric gases, the infrared radiation energy value Lλ measured by the satellite is composed of three main parts: the actual radiative brightness from the ground that passes through the atmosphere to the satellite sensor; upwelling radiance L↑, and descending radiance L↓. The formula is as follows: ε is the land surface emissivity; B (TS) is the blackbody radiance; TS is the LST temperature [K]; and τ is the atmospheric transmittance.
The formula for B (TS) is: The three basic parameters, atmospheric transmittance, upwelling radiance, and descending radiance, must be obtained using the Atmospheric Correction Parameter Calculator in NASA (https://atmcorr.gsfc.nasa.gov/); the imaging time of the remote sensing data and the central latitude and longitude of the downloaded image are input to retrieve the radiance data at that time and place.
The surface emissivity is calculated with the band-math tool in the ENVI software. After estimating the brightness of blackbody radiation corresponding to the surface temperature, the real surface temperature can be obtained according to Planck's law. The formula is as follows: K1 and K2 are two constants whose values can be retrieved from the Landsat 8 files. For band 10 in the TIRS data, K1 = 774.89 W/ (m 2 •µm•sr); K2 = 1321.08 K.

Calculation of normalized vegetation index and specific radiation rate from the surface
This paper adopts the same method for calculation of the surface specific radiation rate as TM/ETM+6 (because the Landsat 8 Band10 band has similar spectral range as the TM/ETM+6 thermal-infrared bands) [19][20]. The specific radiation rate from the surface is estimated based on the weighted mixed model of land-cover types proposed by [ NDVIV is the total vegetation cover index and NDVIS is the vegetation cover index for bare ground. Generally, the simple vegetation calculation model is adopted for calculation.
Empirically, the total vegetation index for our study area is 0.70, and the total bare-ground vegetation index is 0.05. In other words, this expression can be understood as stating that when the NDVI of a pixel in an image is greater than 0.70, Pv is 1, which means that the pixel area is in a region covered by dense vegetation and no bare soil can be seen. When the NDVI is less than 0.05 and Pv is 0, the pixel area is completely bare, without any vegetation coverage.
In the ENVI band-math tool, the formula is: b1 is a numerical image of NDVI, and the image of vegetation coverage (Pv) is obtained after using the band-math tool, with values between 0 and 1. The surface-emissivity image can be obtained with the continuous-input expressions in the band-math tool.

Blackbody radiation and surface temperature calculation
The radiance calculation formula implemented in the band-math tool is: The variables b1 represent the numerical image of the surface's specific emissivity, b2 is the radiation brightness numerical image of the 10th band after calibration, a1 is the brightness of upwelling radiance taken from NASA data, a2 is the transmittance of the atmosphere in the thermal-infrared band, a3 is the descending radiance. The radiance image of a blackbody at the same temperature is calculated.
When all the relevant parameters are obtained, Planck's law can be used to calculate the surface temperature. Since the calculated temperature is in Kelvin temperature, 273.15 is subtracted to return the temperature in Celsius.

Image mosaic
As the study area is not covered by a single Landsat 8 image, multiple remote sensing images must be mosaicked together, and the study area is cropped from the resulting mosaic.
The seamless mosaic scripts built into the ENVI software allow users to combine multiple non-geographic images into a larger single image. The resulting image is not too stiff for radiation calibration and introduces no error to the results.
The simple process of mosaicking is performed by the seamless mosaic tool in ENVI software: 1. Data loading; 2. Color correction; 3. Seamlines/feathering; 4. Export. We especially found that, although the color correction and add-seam-lines operations can make the mosaic images match better and can prevent bad pixels at the contact edges, these two operations will change the original DN values in the images, and this may cause serious errors in the retrieval of surface temperature. Thus, the above operations were not applied.
Taking September 1, 2014, as an example, the result after mosaicking is shown in Figure 6: Obvious contact lines appear at the joints between the images, which are caused by the 1 absences of color correction and seam lines. However, it should be noted that the Landsat 8 2 satellite revisits an area every 16 days, and sometimes the data from two images taken at very 3 close imaging times may be unusable due to too much cloud cover (>50% ). In this case, only 4 images with a relatively long time difference can be used for mosaicking, which leads to 5 obvious contact lines between the images, and discontinuities in the data at these boundaries, 6 but these differences do not seriously affect the final results. 7 As the research area is only a small region of this image, the images are cropped 8 according to the range of the selected research area. We cropped the images according to 9 latitude and longitude.

Verification of LST retrieval results 42
To verify the accuracy of an retrieval algorithm, it is usually necessary to compare the 43 results against values measured on the ground or other standard values. However, due to the 44 difficulty in obtaining these measurements retrospectively and over a large area, the historical 45 13 temperature record of cities and counties provided by the Liaoning meteorological bureau 46 could be used. However, the historical temperature data at the county-level are difficult to obtain and have low credibility. Hence, this paper only used city-level data for analysis. In 48 addition, MODIS LST data were used for verification [21]. 49 Since the spatial resolution of Landsat 8 is 30 meters (the real resolution of the sensor is 50 100 meters) and the spatial resolution of MOD11_L2 is 1km, the Landsat 8 temperature 51 retrieval results need to be translated to the same resolution for comparison [22]. Therefore, 52 Landsat 8 needs to be resampled to 1km resolution to match it to the MOD11_L2 data [23]. 53 Table 1 Table 2: 60 The area of Fengcheng after resampling the Landsat 8 results is shown in Figure.      The data in Table 2 shows that Landsat 8 point positions temperature and the 86 temperature of from the MODIS geothermal dataset were similar, most within ±2ºC of each 87 other, some within ±3℃, and only a few points had a serious difference in temperature. The 88 point of great difference is likely caused by different initial resolutions or the urban 89 heat-island effect and appears in Figure 15, Figure 16, and Table 1. The maximum 90 temperature, the lowest temperature, and the average temperature of the two kinds of data 91