Exploring and monitoring geothermal and volcanic activity using Satellite Thermal Infrared data in TVG, Taiwan

Tatun Volcanic Group (TVG) in the northernmost tip of Taiwan is an active volcano area with the possible magnitude of volcanic explosivity index (VEI) 4 which can devastate the Taipei metropolitan area if eruption occurs. To integrate Thermal Infrared (TIR) remote sensing for exploring geothermal energy and monitoring volcanic activity in TVG, Landsat 7 ETM+ imagery is used to retrieve the Land Surface Temperature (LST) and Radiative Heat Flux (RHF); MODIS 8-day average LST products are adopted for monitoring Jinshan fault zone of TVG. Validated results show that occurrences of hot springs and fumaroles conform to LST anomaly distribution and the overall temperature in E-W ridge is higher than the SW-NE ridge in TVG. Thermal anomaly patterns indicate that distributions of higher LST areas appear correlating with the development of Jinshan fault; for further verification, Hilbert-Huang Transform (HHT) was applied for decomposing MODIS 8-day average LST time series (2002 2016) in Jinshan fault zone. Possible related physical processes underneath Ensemble Empirical Mode Decomposition (EEMD) components based on HHT were also elaborated. The inference that LST component with average period around a month (EEMD component 1, i.e., C1) has irregular spikes which is likely associated with earthquakes in TVG has been investigated in detail. Finally, HHT comparison analysis with three active volcanoes in Philippines and Indonesia is performed for assessing the eruption potential of TVG. Preliminary results show that the regular C4 annual cycle from TVG’s LST time series implies the TVG’s current status is calm and resting with no sign of eruption for future decades. Article history: Received 15 November 2017 Revised 20 January 2018 Accepted 22 January 2018


InTroducTIon
Tatun Volcanic Group (TVG) consists mainly of andesitic lava domes in the northernmost tip of Taiwan.It is within the domain of Yangminshan National Park and the only area with the active volcanic features such as fumaroles and hot springs in the island.Significant geothermal and micro-earthquake activities near the surface occur despite no historical record of eruption.TVG, in between the Jinshan fault and Kanchiao fault (Fig. 1), is built on a late Tertiary sedimentary basement by the multiple eruptions during the Quaternary Ryukyu arc-related volcanism.It includes over 20 well-preserved volcanic edifices built mainly of lava flows and pyroclastic and water-laid volcanoclastic strata.Two eruption periods can be grouped based on the age of volcanic strata: the first period includes small-scale eruptions around 2.8 -2.5 Ma with limited volumes of lava and pyroclastic deposits, whereas the second period occurred around 0.8 -0.2 Ma with the formation of a cauldron, followed by voluminous lava flows at most of the volcanic centers.Mt.Qixing, Mt.Dajianhou, and Mt.Huangzuei are the youngest volcanoes of the TVG with ages of 0.3 -0.2 Ma (Song et al. 2007).A previous study of Helium isotope ratios signifies magma chambers might exist beneath the TVG (Yang et al. 1999).A study from volcanic ashes of TVG indicates that the last eruption with the possible magnitude of VEI 4 (volcanic explosivity index) might be around 5500 year ago (Belousov et al. 2010).Another study also suggests that potential volcanic hazards consists of lava, pyroclastic flows and lahars, if assuming Mt.Qixing as the future eruptive center and erupts in the past manner (Tsai et al. 2010).In addition, seismic observations indicate that a magma chamber may still exist beneath TVG and that a future eruption is possible.Analysis of crustal earthquakes (depth < 30 km) suggests that local seismicity in the TVG seems to be driven by circulation of hydrothermal fluids from subsurface hydrothermal or volcano-related activities.Focal mechanism solutions in TVG are primarily normal faulting, indicating earthquakes are most likely conjoined with hydrothermal and magmatic activities in a back-arc extensional setting (Kim et al. 2005;Konstantinou et al. 2007).Besides, seismogenic stress field beneath the TVG exhibits a localized axes orientation, especially in the Mt.Qixing area, which may be connected with the opening of microcracks and the ascent and circulation of fluids in the upper crust (Konstantinou et al. 2009).
TVG is considered a potentially active volcano based on various observations from previous studies.It may cause the severe social impact in Taipei if there is any potential volcanic activity in TVG.That is, it will devastate the Taipei metropolitan area with approximately 8 million populations if eruption occurs (Konstantinou 2015).Besides, two nuclear power plants (Jinshan and Kuosheng nuclear power plants) located in the proximity (15 km northeast and 12 km east, respectively, of the volcanoes) has also leveled up its significance and imperativeness.
TVG has been designated and reserved as part of the Yangmingshan National Park since 1962.Many volcanic features such as composite and conic volcanoes, fumaroles, craters, and crater lakes are well-preserved in the park.Apparent volcanic activity includes sulfur fumaroles, steam vents, geothermal energy and hot springs are prevalent in the Xiaoyoukeng, Dayoukeng, and Liuhuanggu areas, etc.The TVG comprises E-W and SW-NE volcanic ridges as shown in Fig. 1.Each ridge is around 15 km in length, and has an average elevation of 800 -1000 m above sea level.The SW-NE ridge has much sharper ridge and steeper alpine because of severe erosion.In contrast, the E-W ridge contains much gentler shape and preserve many primary volcanic edifices such as cones, craters, lava flows, and pressure ridges (Yang et al. 2004).Considering the facts that both ridges are located in similar weather conditions, and have similar composition of rocks, these features signify that the E-W ridge is much younger than the SW-NE ridge.The E-W ridge consists of 23 major well preserved volcanic landforms, the highest of which is Mt.Qixing at 1120 m above sea level.Slopes of most volcanic edifices are up to 40 degrees and the size is small-to-moderate.Andesite is the major production of TVG all through its geological history.The SW-NE ridge composes mostly two-pyroxene andesite without hornblende, while the E-W ridge consists of andesite with high rates of hornblende.This implies a huge change in water containing of magma possibly caused by the shifting of volcanic plumbing system and volcanism from the SW-NE ridge to the E-W ridge (Belousov et al. 2010).
Taiwan Volcano Observatory at Tatun (TVO) is the first program for volcano monitoring which was initiated by the Ministry of Science and Technology (MOST) in 2011.Four different approaches were integrated to study the TVG in TVO for the first phase.Namely, soil gas monitoring, crustal deformation survey, magnetotelluric measurement, and geothermal observation (ground stations).Observation results from TVO indicate that volcanic activities in the Tatun volcanic area still active.And temporal variations on micro-earthquakes and geochemical measurements are particularly significant.Magnetotelluric measurements suggest that the possible magma chamber at the shallow crust might be cooling now, but the occurrence of other deep magma sources can't be eliminated.Thus, it is a very urgent task to deliberately investigate and monitor the TVG nowadays.Meanwhile, further multi-discipline studies are encouraged to improve the understanding of the potential activities in the TVG.And the geothermal resource of TVG might be obtained for the civilian purpose in the future (Lin 2009).
Exploring and monitoring the geothermal and volcanic activity is challenging in the aspects of time and cost-efficiency.Integration of TIR remote sensing technology with existing geophysical methods for the volcano monitoring is applicable to facilitate the understanding of TVG.Thus, the application of satellite remote sensing for detecting and quantifying thermal anomalies has been conducted.Satellite-based infrared data both from Landsat 7 and MODIS are utilized in this research.Landsat 7 TIR imagery was used to retrieve the LST and RHF; MODIS LST and Emissivity (LST&E) product was also employed to the HHT data analysis for Jinshan fault zone in TVG from 2002 -2016.Prime objectives of this paper are: (1) to explore LST anomaly of TVG and estimate RHF of TVG, (2) to monitor TVG by applying the HHT for decomposing MODIS LST time series along the Jinshan fault zone in TVG, and (3) to assess the eruption potential of TVG by analyzing corresponding physical processes underneath components of MODIS LST decomposition and comparison to other erupting volcanoes.

MATErIAlS And METhodS
This research uses both Landsat 7 ETM+ LT1 data and MODIS LST products.MODIS LST products are ready for further analyses while Landsat data is required for postprocessing to obtain the LST and RHF results.The nominal spatial resolution for Landsat LST and MODIS LST are of 30 m and 1km, respectively.Landsat 7 ETM+ data are used to retrieve the LST anomaly pattern in TVG, while MODIS LST products are for the monitoring of the volcanoes and Jinshan fault systems.

landsat 7 Products
Landsat program, initiated in 1965 by the US Geological Survey (USGS), was inspired by the Apollo Moon-bound missions with the idea to gather facts about the natural resources of our planet.The first civilian Earth observation satellite Landsat 1 was launched on 23 July 1972.Landsat 7 (launched in 1999) and Landsat 8 (launched in 2013) are currently active missions.Landsat 7 Enhanced Thematic Map-per Plus (ETM+) imagery contains eight spectral bands.The spatial resolution for Bands 1 -5 and 7 is 30 m, Band 6 is 60 m, and for Band 8 (panchromatic) is 15 m.All bands can acquire one of two gain settings (high or low) depending on radiometric sensitivity and dynamic range, while Band 6 (thermal band) collects both high and low gain for all scenes (Arvidson et al. 2001).
The major Landsat scene in this study was obtained on 3 December 2001 and covers northern part of Taiwan (WRS2 Path/Row 117/43; Acquisition time: UTC 02:09:19, local time 10:09:19 AM).The scene size by default is about 170 km north-south by 183 km east-west (106 miles by 114 miles).The image is in good quality of rank 9 (Image Quality ranges from 0 -9 with 9 being the best).Clear weather condition was presented with average wind speed of 2.3 m s -1 and no precipitation according to the meteorological report of Central Weather Bureau (CWB) of Taiwan.As specified in the Landsat 7 Science Data User's Handbook (from National Aeronautics and Space Administration, NASA), Landsat scenes are processed to Standard Terrain Correction (Level 1T -precision and terrain correction) which gives systematic radiometric and geometric accuracy through integrating ground control points and topographic accuracy by the Digital Elevation Model (DEM) (Chander et al. 2009).Spectroradiometer) is the key instrument aboard the Terra and Aqua satellites launched on 18 December 1999 and on 4 May 2002, respectively.Terra MODIS and Aqua MODIS are orbiting the entire Earth's surface to provide global coverage every one to two days, acquiring data in 36 spectral bands.The Aqua MODIS instrument acquires data twice daily (it crosses the equator locally at 1:30 PM and AM), as does Terra MODIS (it crosses the equator locally at 10:30 AM and PM).These four daily MODIS observations serve to advance global monitoring of the Earth.MODIS instruments provides high radiometric sensitivity (12 bit) in 36 spectral bands ranging in wavelength from 0.4 -14.4 μm (VIS, NIR, SWIR/MWIR, and LWIR spectral regions).Two bands are imaged at a nominal spatial resolution of 250 m at nadir, with five bands at 500 m, and the remaining 29 bands at 1 km.A ±55-degree scanning pattern at the NASA's Earth Observing System (EOS) orbit of 705 km achieves a 2330-km swath.

MODIS (Moderate Resolution Imaging
This research used MODIS LST & Emissivity 8-Day L3 Global 1km (Product ID: MYD11A2) nighttime datasets acquired in Bands 31 and 32 of Aqua MODIS from 2002 -2016.The level-3 MODIS global LST and Emissivity 8-day data are composed from the daily 1-km LST product (MY-D11A1) and stored on a 1-km Sinusoidal grid as the average values of clear-sky LSTs during an 8-day period.MYD11A2 is comprised of daytime and nighttime LSTs, quality assessment, observation times, view angles, bits of clear sky days and nights, and emissivity estimated in Bands 31 and 32 from land cover types with the collection of version 5. Version-5 MODIS products are validated to Stage 2, which means that their accuracy has been assessed over a widely distributed set of locations and time periods via several ground-truth and validation efforts (Wan et al. 2004;Wan 2007).

hilbert-huang Transform (hhT)
The combination of the Hilbert spectral analysis (HSA) and the empirical mode decomposition (EMD), designated as the HHT by NASA.The purpose of HHT is to decompose MODIS LST time series into multiple components in the time domain, and each component is related to a specific physical process.The HHT is suitable specifically for analyzing nonlinear and nonstationary data.The key part of HHT is EMD with which any complicated data set can be decomposed into multiple intrinsic mode functions (IMFs).Each IMF represents a narrow band frequency-amplitude modulation that is often related to a specific physical process.This decomposition method is adaptive, data-driven and highly efficient.It is applicable to nonlinear and nonstationary processes as the decomposition is based on the local characteristics of the data (Huang and Wu 2008;Wang et al. 2014).
EMD is the empirical and adaptive algorithm for the accurate expression of a time series in the time-frequencyenergy domain.In EMD, the data x(t) are decomposed in terms of IMFs, c j : where and r n is the data residual of x(t) after IMFs have been extracted n times.The EMD is implemented through a sifting process that uses only local extrema.For any data set, x(t) = r j -1 , the procedure is as the following 3 steps: (1) take all the local extrema (both maxima and minima) and connect them with a cubic spline as the upper (lower) envelope and calculate the corresponding local mean; (2) extract the first component h by subtracting local mean from the data x(t); and (3) treat h as the new data and repeat steps 1 and 2 until the envelopes are symmetric about zero to a certain threshold.The final h is designated as c j .When the residue r n becomes a monotonic function or no more IMFs can be extracted from the above procedure, the sifting process is finished.The secular trend of a time series is then obtained after all the oscillatory components (riding waves) are removed by applying EMD.EMD has accumulated thousands of citations by the numerous successful applications in vari-ous scientific and engineering fields ever since its development about 17 years ago (Wu et al. 2011).However, one of the major drawbacks of the original EMD is the frequent appearance of mode mixing.To overcome this shortcoming, a noise-assisted Ensemble EMD (EEMD) method is proposed.EEMD is a substantial improvement over the original EMD.EEMD includes sifting an ensemble of white noise-added signal and deals with the mean as the final result.As the EMD is a time domain analysis method, the white noise is averaged out with sufficient number of trials; the only remaining part survives the averaging process is the signal, which is then served as the true and more physical meaningful answer (Wu and Huang 2009).

data Processing for landsat 7 Products
LST retrieval is based on the Blackbody Radiation and the Planck Function.Planck Function is used to calculate the radiance emitted from a Black Body.Brightness temperature is derived by the inverse of the Planck Function.Landsat 7 ETM+ sensor acquires radiance information and stores it to the Digital Number (DN) between 0 and 255 (8bit data).So the LST can be retrieved by converting these DN values to degrees Kelvin.In this study, main steps of data processing include radiometric calibration, atmospheric correction, topographic correction, and emissivity calculation.Specific procedures are illustrated as shown in Fig. 2 (Sobrino et al. 2008;Chan et al. 2018).
Radiometric calibration is the first step to convert the DN value into the Top of Atmosphere (TOA) radiance according to formulas from the Landsat 7 Science Data Users Handbook.Two atmospheric correction modules are adopted to Landsat scene depending on the wavelengths.For wavelengths in the visible through near-infrared and shortwave infrared regions (up to 3 μm), Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) is applied.For the thermal band, a web-based atmospheric correction tool "The Atmospheric Correction Parameter Calculator" provided by NASA for single thermal band sensors is used; the site-specific atmospheric transmission, and upwelling and downwelling radiances can be calculated from the calculator (Barsi et al. 2003).
Emissivity calculation is based on the work of Sobrino et al. and referred to as NDVI threshold methods (Sobrino 2001).It uses the general method to incorporate emissivity as a function of the atmospherically corrected red and nearinfrared bands.
The LST result is computed from a single-channel algorithm proposed by Artis and Carnahan (1982).Based on LST, theoretical RHF can be calculated using the following Stefan-Boltzmann equation (Bromley et al. 2011;Mia et al. 2012), where Qr = radiative heat flux (W m -2 ), x = atmospheric transmissivity, v = Stefan-Boltzmann constant, f = emissivity, T s = land surface temperature (K), and T a = ambient temperature (K).Stefan-Boltzmann law states principally the total radiating energy of a blackbody, which is immersed in a medium with ambient temperature (i.e., absolute temperature of the surrounding medium).For the RHF calculation based on LST, the near-surface air temperature serves as the ambient temperature.However, it is not feasible to acquire the two dimensional air temperature distribution right at the particular date and time of Landsat 7 satellite's overpass, thus the air temperature from the Central Weather Bureau (CWB) meteorological Zhuzhihu station's hourly data was interpolated (to Landsat overpass date and local time 10:09 AM) and used as the ambient temperature.

lST and rhF distribution of TVG from landsat 7 ETM+
Figure 3 shows the retrieved LST distribution of TVG from Landsat 7 ETM+ band 6 (High Gain setting, band 6H or 62) data on 3 December 2001.The pattern indicates that the lowest LST in the study area is about 18°C and the highest is 34 °C with color spanned from blue to red.The squared area shows the thermal anomaly at E-W and SW-NE volcanic ridges of TVG.High LST areas are overall 3 -6°C higher compared to the ambient background temperature.The pattern of thermal anomaly is consistent with the occurrence of hot springs, which are indicated by the black dots.Results imply that occurrences of hot springs and fumaroles area are in agreement with anomaly areas.The overall LST in E-W ridge is generally hotter than the SW-NE ridge, for mostly hot springs and fumaroles area are located at E-W ridge.Furthermore, thermal anomaly patterns also indicate that distributions of higher LST areas appear correlating with the development of Jinshan fault in TVG as shown in the polygonal Region of Interest (ROI).Figure 4 illustrates the RHF of squared area (which is less contaminated by the anthropogenic activity) in TVG based on LST from Landsat 7 ETM+ according to Eq. ( 2).The RHF have the similar anomaly pattern with LST.And the max and mean value of RHF are 52.7 and 11.9 W m -2 , respectively.

Topographic correction of lST and lST comparison at hot Spring Sites and the Proximal 1 km radius Area
For minimizing the LST changes caused by the rugged terrain in TVG, topographic correction is performed to remove the topographic effects.The empirical Lambertian cosine correction (C-Correction) based on the 20-m-resolution DEM (from the Dept. of Land Administration, Ministry of the Interior, Taiwan) is applied for the topographic correction in this study.The Lambertian model defines that reflected radiance from the inclined surface (L T ) is related to the horizontal surface (L H ) by the following equation: where i is the incidence angle (i.e., the angle between the illumination and the normal to the surface) and z is the sun zenith angle (i.e., the angle between the vertical and the sun).
Topography corrected pattern of LST distribution in TVG are shown in Fig. 5.The positive anomalies LST is well correlated with the location of hot springs in volcanic range which is rural area.In the contrast, hot spring sites located in or near the urban area are contaminated by the anthropogenic activity (i.e., LST anomalies in Taipei City and New Taipei City).Generally, hot springs are characterized by higher LST.To verify the observation quantitatively, zonal analysis in the circular 1-km radius proximal zones centered at hot spring sites were conducted.The comparison of LST at hot spring sites and the surrounding area are as shown in Table 1.

cross Validation of landsat derived lST with ModIS lST Product
Landsat 7 ETM+ dataset has the advantage of high spatial resolution (60 m, resampled to 30-m pixels) but the low temporal resolution (revisit in 16 days) is a disadvantage.Compared to Landsat 7, Terra MODIS provides daily imagery with spatial resolution of 1 km in the study area.According to the general accuracy statement from MODIS land team.The LST accuracy is better than 1 K (0.5 K in most cases), as expected pre-launch (Coll et al. 2009).Thus MODIS LST product is used for the cross validation of Landsat retrieved LSTs.MODIS/Terra LST Daily L3 Global 1 km (Product ID: MOD11A1) daytime LSTs acquired by Terra MODIS at local time 10:30 AM on 3 December 2001 as displayed in Fig. 6.The comparison on the statistics of retrieved Landsat LSTs and MODIS LST products in TVG is shown in Table 2.The differences on LST of two sensors shown in Table 2 can be caused by various factors, such as the MODIS 21-min delay retrieval after Landsat ETM+.The temperatures may rise rapidly as the time approaches toward the noon and the highest temperatures are present at the noon.The scale-mismatch effect between the 30 m and 1 km spatial resolutions is also considerably responsible for discrepancies (Chan et al. 2018).

ModIS 8-day Average lST Products (2002 -2016)
for Monitoring TVG Landsat 7 dataset is not sufficient for long term monitoring of the LST in TVG due to the satellite's low temporal resolution (revisit in 16 days).The imagery is often covered    by the cloud.However, the high spatial resolution (60 m) is the advantage.Compared to Landsat 7, MODIS provides daily imagery with spatial resolution of 1 km in the study area.So the MODIS dataset since 2002 is used for the time series analysis in TVG, especially at the Jinshan fault.Jinshan fault is the primary reverse fault (however, it shows extensional tectonism in TVG) caused by the tectonic stress during the orogeny process, and has been wildly studied in previous studies (Chu et al. 1998;Liu et al. 2007).MODIS nighttime data is used to minimize the solar heating effect for the analysis.Figure 7 shows the MODIS 8-day average LST at Jinshan fault zone of TVG (specified polygonal ROI in Fig. 3) from 2002 -2016.MODIS LST time series includes rich characteristics with physics behind the data.To extract the component of different periods of LST.HHT is used to decompose the data.The MODIS LST data and its EEMD are displayed as in Fig. 8.The top plot is the original MODIS LST; C1 -C8 are the decomposed components, and the bottom plot is the trend of LST.
Each EEMD component is often related to a specific physical process.The goal is to search for the corresponding physical process underneath available.Obviously, The C4 with period 370 days (around one year) is in good agreement with the dominant frequency of LST, which represents the annual variation of LST caused by the Sun, season and vegetation, etc. (Chen et al. 2006).C1 and C2 has similar amplitude with average period around a month and two months, respectively.C1 has irregular spikes superimposed on random high-frequency oscillations.These spikes may be related to earthquakes in TVG (Ma et al. 2006;Chen et al. 2013).C3 has an average period one-third year, with a relatively smaller averaged amplitude may be caused by the atmospheric circulation and other factors (Huang and Wu 2008).C5 -C8 have averaged amplitude 1 order smaller than other components.C5 and C6 are the variations of interannual timescales.C5 is quasi-biannual, and C6 has an average period around half a decade.C7 and C8 occurring in a decadal times scale.C7's period is one and half decade, and C8 has period of few decades and even longer.The trend of LST in the bottom plot shows very little variation between 2002 and 2016 but still in an upward trend, which may indicate the volcanoes in TVG are relatively calm.90°C and those in Dayoukeng can even reach 120°C.The overall LST in E-W ridge is generally higher than the SW-NE ridge.This feature verified the proposal from previous studies (Chen et al. 2003;Yang et al. 2004;Liu et al. 2007;Belousov et al. 2010) that the E-W ridge is much younger than the SW-NE ridge.For the younger E-W ridge conserved most volcanic landforms while the older SW-NE ridge preserved none.Thermal anomaly patterns indicate that distributions of higher LST areas appear correlating with the development of Jinshan fault in TVG.To verify this observation, Landsat multi-temporal Brightness Temperature imagery (6 scenes) for the verification of the LST anomaly results is performed as amplified in next section 4.2.And to justify this argument, HHT was adopted for analyzing MODIS 8-day average LST time series in Jinshan fault zone of TVG as elaborated in section 3.4 and 4.3.

landsat Multi-Temporal brightness Temperature Imagery for the Verification of the lST Anomaly results
LST anomaly is affected by soil moisture in different seasons of the year based on the empirical assumption that soil moisture is the main source of variation for LST (Sandholt et al. 2002).However, provision of soil moisture data for analysis is particularly problematic in this study.A single Landsat ETM+ scene on 3 December 2001 maybe not sufficient for the proof of LST anomaly results shown in Fig. 3.For the verification of LST anomaly in TVG.Addi-tional 6 imageries of different year and season (Acquisition date: 2003(Acquisition date: -11-07, 2004(Acquisition date: -11-09, 2008(Acquisition date: -04-26, 2009(Acquisition date: -05-15, 2012(Acquisition date: -05-07, and 2015-06-17) -06-17) which selected meticulously from the Landsat archive were processed to obtain the Brightness Temperature (BT) of the study area as shown in Fig. 9 and Table 3.The BT is the radiance measurement of the Electromagnetic radiation traveling upward from the top of the atmosphere to the satellite, expressed in units of the temperature of an equivalent black body.The surface emissivity is the ratio of the measured spectral radiance over the theoretical blackbody spectral radiance for a surface.It can be converted to the ratio between the BT and LST.Since surface emissivity in the nature is smaller than one.So the LST is always higher than BT.In the perspective of data processing, BT retrieving does not require complicated scene-specific atmospheric parameters and emissivity calculation (as shown in the flowchart of Fig. 2) but provides the similar temperature pattern as LST in the local scale.
All the selected images are suffered from the cloud coverage or image gap caused by the satellite mechanical malfunction.However, the study-interested TVG ridges are mostly retained in the imagery.The imagery set from 2003 -2015 shows clear consistency in the thermal anomaly area discussed above (as shown in Figs. 3, 4, and 5.) which provides ample support for the LST anomaly results.

ModIS 8-day Average lST Time Series
C1 has irregular spikes can be related to volcanic   earthquakes in TVG.Since volcanic earthquakes are caused by the movement of magma.The moving magma exerts pressure on the rocks to the point to crack it.Then the magma gets into the crack and rebuilding pressure.Every time the cracking rock causes a small earthquake.These earthquakes are generally too weak to be felt by human but can be detected and recorded by sensitive seismographs such as in TVG.
Heat transfer in geothermal systems realized generally from geothermal reservoirs (i.e., magma chamber) which have numerous near-vertical faults and relatively impermeable intrusive interspersed in the aquifers.And the thermal anomalies can be detected by surface manifestations, aerial or spaceborne infrared surveys, geochemical analyses, or exploratory drillings (Cheng 1979;Zhao et al. 2008).Primary extensional tectonism in TVG was shown by proofs of seismology and structural geology.A NE-SW narrow fault zone (in a size about 3.5 km width by 16 km length) parallel to the surface trace of Jinshan fault in TVG was delineated from the digital terrain model and field validation.The zone encloses both tectonic fractures and hydrothermal activities in the TVG (Chu et al. 1998).In addition, the possible modal of the volcano-hydrothermal system beneath the Mt.Qixing-Dayoukeng area has been constructed based on the local geology and the seismicity characteristics (Konstantinou et al. 2007).In this regards LST of the Jinshan fault zones in TVG can be the viable heat release channel of the magma chamber beneath and serves as the observation window for the volcanic activities.
For detail investigation of this inference.Earthquake catalog (Epicenters distributed within 15 km centered at Mt. Qixing) from Central Weather Bureau (CWB) of Taiwan during 2002 and 2016 is adopted for the comparison analysis.The earthquake magnitude is related to the energy radiated from the earthquake source.Magnitude and energy has a logarithmic relationship and energy can be derived from the Gutenberg-Richter magnitude-energy formula (Gutenberg 1956).Volcanic activities (i.e., movement of magma) underneath TVG can be the major cause of earthquakes.Seismic evidence has approved the existence of a magma reservoir beneath TVG in the shallow crust (Lin 2016).Thus, the released energy from earthquakes can be responded by LST along the fault area in TVG.
Figure 10 shows the local magnitude-depth plot and magnitude-time plot in TVG from CWB earthquake catalog.Most events happened at a shallow (crustal) depth less than 30 km and are primarily volcanic in origin.Earthquakes detected by the CWB increase suddenly in numbers after 2012.The reason is that CWB replaced and updated all 16-bit (resolution) seismographs to 24-bits in 2012.To focus on the volcanic origin events and recede the effect of detectability inconsistency after 2012, earthquakes with depth more than 30 km or magnitude less than 1.5 are discarded for further analysis.For facilitating the comparison, magnitude-time dataset of earthquake catalog was first re-sampled into 8-day summing magnitude to align with the LST dataset.As the result, Comparison between the HHT decomposition of LST and 8-day summing magnitude-time plot in TVG is displayed in Fig. 11, which is an enlargement of C1 of Fig. 8.
For the comparison of LST C1 and 8-day summing magnitude datasets.Cross correlation to measure the similarity of two series was applied as shown in Fig. 12. Referring the perfect correlation of two time-series shows a single spike at the zero lag in the cross-correlation plot.It appears fair correlation between the LST and earthquake 8-day summing magnitude series, which signifies the possible correlation between the C1's spikes and earthquakes in TVG.Heating from circulation of hydrothermal fluids during the earthquakes may cause the increase of LST along the fault structure in TVG.Since seismicity in the TVG was mainly driven by circulation of hydrothermal fluids from the magma chamber beneath.In addition, previous study of thermal satellite data application in various areas around the globe shows that thermal anomalies appeared about 6 -24 days before and continued about a week after earthquakes (Tronin 2006).Delay time of the top three correlation values shown in Fig. 12 are 35, -46, and -36 days, respectively.Delay of -36 days means LST anomaly appeared 36 days before earthquake occurring, which is loosely agreeable with previous global results.

hhT comparison Analysis with Active Volcanoes in Philippines and Indonesia
To assess the eruption potential of TVG, true physical processes underneath the EEMD components (C1 -C8) are quests to reveal.However, they are generally hard to obtain owning to the unknown natural complexity of volcanic systems.The feasible way to understand the eruption potential of TVG is comparing with the current active volcanoes around the world.Considering both Philippines and Taiwan locate in the tectonic boundary between the Eurasian Plate and the Philippine Sea Plate, Mayon volcano (Stratovolcano; summit elevation 2462 m) and Canlaon volcano (Stratovolcano; summit elevation 2435 m) of Philippines are selected for this end, plus Paluweh volcano (Stratovolcano; summit elevation 875 m) of Indonesia.The MODIS night time 8-day average LST time series for these 3 active volcanoes were retrieved from craters and decomposed by HHT as shown in Figs. 13,14,and 15,respectively.HHT analysis of LST at TVG in Fig. 8 compared with Figs. 13, 14, and 15, a few observations on the LST and EEMD components can be noticeable: (1) the decomposition C4 (Component 4 with annual period) of LST from currently active volcanoes' obviously deviated from the regular annual cycle.This phenomenon is caused by the heat of volcanic eruptions.It implies that the deviation of C4 from regular annual cycle can be an indicator of the level of restless of      active volcanoes.Thus, the fact that there is a regular C4 annual cycle from TVG's LST time series indicates the TVG's current status is quite calm and resting.(2) C1 and C2 are responsive to the volcano eruptions, this may be due to the nature of short period of these components (one to two month).
(3) LST appears an increasing tendency while the eruption occurs in Mayon volcano and Paluweh volcano.However, the match between the eruptions and high LST is not consistent in all cases.One reason is the missing of dataset during the period from 2002 -2016.In the case of Canlaon volcano.For example, there should be 679 data points for the full dataset.However, due to the weather condition or other reasons, 249 data points are missing.Sometimes, the missing data occurred just before or after the eruption events, sharpness of the curve usually shows sparse observation data points in the LST time series.

concluSIonS
This study aims to apply and integrate TIR remote sensing technology for exploring geothermal energy and monitoring volcanic activity in TVG of Taiwan.LST and RHF have been retrieved from the calibrated NASA Landsat 7 ETM+ data.Geothermal anomaly areas are validated by the field investigation of fumaroles, steam vents, and hot springs.MODIS daily LST product is used for the cross validation of Landsat derived LST.Landsat multi-temporal brightness temperature imagery (6 scenes) for the verification of the LST anomaly results were also conducted.Results imply that occurrences of hot springs and fumaroles areas are in agreement with retrieved LST anomaly and RHF areas.And thermal anomaly pattern indicates that distributions of higher LST areas appear correlating with the development of Jinshan fault which is a narrow fault zone with extensional tectonism.
To validate this argument, HHT was adopted for analyzing MODIS 8-day average LST time series in Jinshan fault zone of TVG.Possible related physical processes underneath of EEMD components (C1 -C8) were discussed.Namely, the C4 with period 370 days (around one year) is in good agreement with the dominant frequency of LST, which represents the annual variation of LST caused by the Sun, season and vegetation, etc.It appears fair correlation between the LST C1 and earthquake 8-day summing magnitude series, which signifies the possible correlation between the C1's irregular spikes and earthquakes in TVG.The LST spikes along the Jinshan fault zone in TVG are most likely caused by the heating from hydrothermal fluids during the earthquakes associated with circulation of hydrothermal fluids from the magma chamber beneath.Finally, based on the HHT comparison analysis with three active volcanoes in Philippines and Indonesia, it shows that the deviation of C4 from regular annual cycle can be an indicator for the restless level of active volcanoes.Consequently, the regular C4 an-nual cycle from TVG's LST time series indicates the TVG's current status is quite calm and resting.
The effectiveness of satellite TIR remote sensing in retrieval of LST using Landsat images was presented.The usefulness of HHT data analysis method on MODIS LST time series was also demonstrated.This work suggests that TIR remote sensing is an asset for mapping and quantifying surface features of the fumaroles area as well as volcanic monitoring.TIR remote sensing provides a rapid way of mapping and quantifying surface features to facilitate the exploring volcanic activity.Compared with conventional in-situ methods of point measurements, remote sensing methods are able to produce reliable results which are spatially representative at larger regional scales.However, it is limited to detecting the surficial and the shallow buried geothermal features.Thus integration of geology, geophysics, geochemistry and remote sensing will be the key to better understanding TVG in the future.

Fig. 4 .
Fig. 4. Illustration of the RHF of TVG based on LST from Landsat 7 ETM+ imagery.Note that the high RHF appeared at the southwest corner might be contributed by the anthropogenic activity of Taipei City.

Fig. 5 .
Fig. 5. Topography corrected pattern of LST distribution in TVG from Landsat 7 ETM+ on 3 December 2001.Circular red 1-km radius proximal zones centered at hot spring sites are displayed for the statistical analysis shown inTable 1.
Fig. 7. Illustration of the LST from MODIS at Jinshan fault zone of TVG from 2002 -2016.

Fig. 8 .
Fig. 8. EEMD decomposition of the MODIS LST dataset from 2002 -2016.The top plot is the MODIS LST data, C1 -C8 are the decomposed components, and the bottom plot is the trend.

Fig. 9 .
Fig. 9. Pattern of Brightness Temperature distribution in TVG from Landsat 7 ETM+ on 7 November 2003, 9 November 2004, 26 April 2008, 15 May 2009, 7 May 2012, and 17 June 2015, respectively.Blank areas (white color) indicate "No Data" caused by either cloud coverage or image gaps.Thermal anomalous areas are illustrated by red color in the imagery.Extent of the square area is less contaminated by the anthropogenic activity has been outlined in Fig. 3.

Fig. 10 .
Fig. 10.Magnitude-depth plot (top) and magnitude-time plot (bottom) in TVG from Central Weather Bureau (CWB) earthquake catalog.Earthquakes with epicenter distributed within 15 km centered at Mt. Qixing were selected and most of the magnitude are smaller than 4. Earthquakes with magnitude less than 1.5 or depth more than 30 km are omitted in this figure and discarded for further analysis.

Fig. 12 .
Fig. 12. Cross correlation of the first EEMD component C1 and earthquake 8-day summing magnitude series.The value of vertical axis in the correlogram indicates the similarity and time delay of the peak value is 35 days.It means delay LST by 35 days to match with earthquake occurrence.Time delay of the secondary and tertiary high values (almost equivalent to the peak value) are -46 and -36 days, respectively.

Fig. 13 .
Fig. 13.EEMD decomposition of the MODIS LST dataset from 2002 -2016 for Mayon crater of Philippines.The top plot is the MODIS LST data, C1 -C7 are the decomposed components, and the bottom plot is the trend.Vertical dashed red lines indicate the recorded eruptions and volcanic earthquake events.

Fig. 14 .
Fig. 14.EEMD decomposition of the MODIS LST dataset from 2002 -2016 for Canlaon crater of Philippines.The top plot is the MODIS LST data, C1 -C7 are the decomposed components, and the bottom plot is the trend.Vertical dashed red lines indicate the recorded eruptions and volcanic earthquake events.

Fig. 15 .
Fig. 15.EEMD decomposition of the MODIS LST dataset from 2002 -2016 for Paluweh crater of Indonesia.The top plot is the MODIS LST data, C1 -C8 are the decomposed components, and the bottom plot is the trend.Vertical dashed red lines indicate the recorded eruptions and volcanic earthquake events.

Table 1 .
Statistics of circular 1-km radius proximal zones centered at hot spring sites.On-Site column refer to the estimated LST right at the site spot.
Note: * Hot spring sites located in the urban area which tend to be contaminated by anthropogenic activity.

Table 2 .
Comparison on the statistics of retrieved Landsat LST and MODIS LST products in TVG.The feature-based (pixel statistics) comparison rather than pixel-by-pixel are adopted for the cross-validation of two LST imagery difference.

Table 3 .
Statistic summary of Landsat ETM+ multi-temporal brightness temperature imagery in the study area.