Thermal pattern at kueishantao (KST) volcano of Taiwan from satellite-observed temperatures and its implication

Kueishantao (KST) is an active volcanic island off the northeastern coast of Taiwan. Tectonically, it lies in the south of the Okinawa Trough and opposite to the Ilan plain, in which is the southwestern end of the trough. KST provides a convenient observation site for the subsurface geological and geothermal activity and mechanism at its proximity. Land surface temperature (LST) of volcanoes detected from satellite sensors reflects the thermal status of heat sources in the subsurface. LST thus is a key parameter to the understanding of the volcanic process and geothermal resources. This research utilizes the satellite-observed multi-temporal land surface temperature imagery in 1999–2022 on the Kueishantao volcano of Taiwan to explore its geothermal state. The U.S. NASA Earth-observing satellites onboard three thermal sensors (i.e., Landsat ETM+, Terra ASTER, and Aqua/Terra MODIS) derived time series of land surface temperature from 1999 are employed to define the past and current pattern of geothermal activity plus the future trend of the KST. The spatiotemporal LST distribution of KST volcano is explored and analyzed. The spatial LST distribution of the KST volcano indicates that LST anomaly areas are mainly located on the southeast island, which is well correlated with the possible magma reservoir location from previous geophysical and geological surveys. An increasing trend of two-decade LST time series is revealed from all three thermal sensors. The retrieved surface thermal pattern shows non-linear temperature variations that imply the non-steady-state nature of the subsurface thermal sources at this volcano. In summary, satellite LST observations facilitate the understanding on the subsurface magmatic processes of active volcanoes for further management of geothermal resources.


Introduction
Kueishantao (KST) is an active volcano with the last volcanic eruption occurring around 7 ka in the Holocene Epoch (Chen et al 2001). It consists of andesitic lava flows and pyroclastic rocks with a summit of 401 m above sea level. It is also known as Guieshan Island, i.e., Turtle Mountain Island for its turtle-like shape with the head facing the Pacific Ocean and the tail pointing to the Taiwan mainland (figure 1). KST is located offshore of NE Taiwan, a 12-km distance from the main island of Taiwan, and symbolizes one of Taiwan's two primary volcanic centers (the other being the Tatun Volcano Group). KST island lies in the southern Okinawa Trough and opposite the Ilan plain, which is the southwestern end of the trough and is formed by back-arc spreading with active microseismicity plus volcanism (Lee et al 1980, Letouzey and Kimura 1986, Sibuet et al 1998, Chung et al 2000. Subsequently, KST provides a convenient observation site for the subsurface geological activity and process around its proximity (Konstantinou et al 2013). In addition, geothermal energy is transferred to the Earth's surface via earthquake or volcanism because both volcanoes and earthquakes have resulted from the Earth's tectonic movement. Understanding the thermal condition and dynamics of a volcanic system is crucial to comprehending the reactivation and persistence of volcanic eruptions (Wadge 1982, Slezin 2003, Mann and Stix 2010, Rosi et al 2022. Land surface temperature (LST) isa key factor in the prediction of eruption and monitoring of an active volcano. The ground's surface temperatures do reflect the temperature underneath it. Thermal features on volcanoes, such as steam vents, geysers, hot springs, lava flows, or lava domes, generally experience fluctuations in surface temperature before a volcanic eruption. It can be useful to foresee changes in activity by identifying these 'thermal anomalies'. This research intends to explore the KST volcano's land surface temperature characteristics and its subsurface status quo by using satellite observations.
Historical records in China's Qing Dynasty (1775-1795 AD) describe the KST island erupting with the red lava flow (GVP 2013). Numerous submarine volcanoes located along the southeast of KST indicate the presence of hydrothermal vent systems (Lee et al 1980, Sibuet et al 1998, Li et al 2016, Lebrato et al 2019. The active submarine fumaroles and sulfatases caused the discoloration of seawater surrounding to make the water surface colored milky white, also known as the 'Milky Sea' (Chen et   eruption likely happened around 7 ka, according to a reported thermoluminescene (TL) age for the bottom portion of the volcanic sequence. The bubble gas samples taken from the sea surface and submarine at KST reveal significantly high 3He/4He ratios (7.35 to 8.39), indicating that the area's magmatic activity is active (Chiu et al 2010). In addition, the earthquakes in proximity have sparked the possible eruption of the volcano since 2014, for the increased seismicity and shallow earthquakes nearby have caused landslides and even hot-mud spiting on this island (CWB 2022). The potential and threat for erupting of the KST volcano are thus everpresent and consistent monitoring for volcanic activities at KST is demanded.
Risk mitigation and resource management require constant monitoring and evaluation of the volcanic activity at KST. However, because such studies demand long-term volcanic observation, quantitative analyses of long-period eruptive behavior are often limiting. Fortunately, remote sensing offers useful data for revealing potentially hidden volcanic mechanisms. We make use of a massive satellite imagery dataset at KST from the Landsat satellite and Terra/Aqua satellites to better understand the past behavior and future eruption trends of this volcano. Specifically, we are conducting spatio-temporal monitoring analyzes from over 20-year time series of land surface temperature obtained from Terra/Aqua and Landsat thermal infrared imagery archive to investigate mechanisms of this volcano and to facilitate effective mitigation of potential hazards and management of geothermal resources.
Satellite thermal remote sensing is a prospective method of monitoring active volcanoes around the world. Compared to in situ observing and on-site monitoring, satellite remote sensing offers greater advantages like the capacity to gather data over vast areas and to monitor surface objects regularly over time. There are currently dozens of Earth-observing satellites orbiting around the Earth. The increasing number of on-orbit satellites subsequently collect larger and larger amounts of observation datasets. For instance, NASA's Earth science data archive (by 2021) is about 40 petabytes in total (1 petabyte is about 1,000 terabytes; 1 terabyte is around 1,000 gigabytes-the storage of 250 feature-length films) (NASA 2022b, NASA 2022c). Given that a large amount of freely accessible satellite thermal datasets (e.g. Landsat, MODIS, and ASTER) are available, the application community and imagery users are expected to move from the stage of image validation to the stage of diverse, sophisticated, and integrated applications (Stehman and Foody 2019, Wulder et al 2019, Bauer 2020. In light of this, we aim to utilize a multi-spatio-temporal dataset to investigate the study area to better comprehend the interactions of geothermal energy between the subsurface activities and the land surface temperatures.

Materials and methods
A thermal remote sensing system generally uses a detector (sensor) to collect the thermal infrared energy that the Earth's surface emits. Sensors record the energy as voltage, which an analog to digital converter turns into the scaled integers (called the Digital Number, or DN) for the energy (Harris 2013). In a similar way, NASA satellites onboard thermal sensors acquire radiation energy from the land surface and converts the energy into the DN on the imagery, i.e., the grayscale images. So the temperature can be converted from DN values to degrees Kelvin by the inverse of the Planck Function (i.e., Planck's radiation law).
We have collected satellite thermal infrared images on the KST volcano for the period 1999-2022 from the NASA imagery archive in this study. Landsat images (30 m spatial resolution), ASTER images (90 m), and MODIS land surface temperature products (1 km) were selected. The selection criteria of these three datasets are based on the spatial and temporal resolution of imagery on KST. In practice, imagery of the single satellite sensor may be inadequate for volcanic monitoring purposes because of its spatial and temporal limits. Thus we have utilized the imagery of three thermal sensors. Fortunately, a complementary relationship exists among Landsat/ ASTER and MODIS imagery-Landsat/ASTER has a high spatial resolution (30 m; 90 m) but a low temporal resolution (16 days); MODIS has a high temporal resolution (1 day) but a low spatial resolution (1 km) (Bindhu et al 2013, Weng et al 2014, Bonafoni 2016).

Landsat imagery retrieved land surface temperature
The thermal infrared bands on Landsat's sensors are in the 10.40-12.50 μm wavelength range. These wavelengths, which are utilized only to detect long-wavelength radiation emitted from the Earth and whose intensity mostly relies on surface temperature, are much beyond the range of human eyesight. The sensors' construction and operation are based on NASA-developed quantum mechanics principles (NASA 2022a). Imagery from Landsat satellites and sensors, specifically Landsat 7 ETM+, are used in this study. The Landsat satellites pass overhead at around 10.30 a.m. local time. The sensor's spatial and temporal resolutions are 30 m and 16 days, respectively (NASA 2019).
To determine the LST on KST, we have chosen thermal imagery from different years and seasons depending on image quality (mostly cloud cover). Figure  have provided more thorough explanations of the Landsat data processing processes and result validation. The Planck Function, which links the rate at which a surface radiates energy to a function of kinetic temperature, is the basis for the single-channel algorithm used to calculate LST (Artis and Carnahan 1982). LST retrieval accuracy is typically better than 2 K (Vidal et al 1997, Barsi et al 2003. 2.2. MODIS 8-day daytime and nighttime temperature imagery dataset MODIS imagery has a shorter revisit period (sub-daily) than Landsat imagery, albeit with a coarse (1 km) spatial resolution. The MODIS sensor is equipped on NASA's Terra satellite and Aqua satellite. Terra MODIS and Aqua MODIS are launched in 1999 and 2002, respectively. Both Terra MODIS (overpass time 10:30 pm) and Aqua MODIS (overpass time 1:30 am) provide the images in the study area. We have incorporated both Terra and Aqua to improve the amount of data. Time series analysis was performed using data from the MODIS 8-day day/night temperature data (Product ID: MOD11A2; MYD11A2) for the years 1999-2022. The MODIS thermal infrared bands (wavelength: 3.6-12.3 μm) are used to produce the 8-day average LST of the corresponding daily LST pixels, which is provided by the MOD11A2/MYD11A2 temperature products (Wan andLi 1997, Wan 2014). According to the NASA MODIS land team, the LST product's accuracy is typically better than 1 K (0.5 K in most circumstances). The exemplary grayscale MODIS, ASTER, and Landsat LST imagery in the KST volcano are shown in figure 3.

ASTER temperature imagery dataset
The ASTER is the thermal sensor onboard NSAS's Terra satellite launched in 1999. It provides land surface temperatures with a temporal resolution of 16 days and a spatial resolution of 90 m. We have used ASTER L2 (level 2) surface kinetic temperature data for cross-validation because it has similar spatial and temporal resolutions with the Landsat imagery. The five Thermal Infrared (TIR) bands in the 8-12 m spectral range are used to extract the Terra/ASTER surface kinetic temperature (Product ID: AST 08; for the land areas only). Land

Ensemble empirical mode decomposition (EEMD)
This work utilizes the EEMD to decompose the LST time series into several simple components, such as the simple oscillatory mode, in the time domain. The original LST time series are fully and orthogonally represented by the resultant LST decomposition components (Huang et al 1998, Huang and Wu 2008, Nunes and Deléchelle 2009, Wang et al 2014. Thus, the trends of the LST time series are constructed by superimposing several simple EEMD components. The EEMD is an effective empirical signal processing technique to handle nonstationary and nonlinear signals and extract useful information from them. (Huang et al 1998, Huang et al 2003). Most traditional signal processing techniques are based on linearity and steady-state assumptions. However, dealing with nonstationary and nonlinear signals at the same time is often needed in real-world circumstances. Specifically, for the alternative relations in the time, frequency, and energy domain, the EEMD is capable of handling nonstationary and nonlinear signals. EEMD is the adaptive and flexible data decomposition method. Its analysis basis is based on the nature of the data and derived from the data. In contrast, many traditional methods of data decomposition are not adaptive (Huang and Wu 2008); they have a presumptive base (i.e., trigonometric functions in Fourier analysis and base wavelet in Wavelet transform, etc). In practice, we use EEMD to extract physically meaningful information from data and to better explain physical phenomena and resolve engineering problems

Results
3.1. Monitoring KST volcano with the 8-day average MODIS LST imagery dataset (1999 to 2022) MODIS 8-day average daytime plus nighttime LST time series at KST has been used for the analysis. The total amount of MODIS observation dataset used in this study is around 3 gigabytes. The original retrieved MODIS LST time series contains large negative values (such as −273.14999°C; the absolute zero temperature) that are caused by the meteorological clouds and other factors. So, we have to screen out data points that are cloud contaminations marking them as NaNs. There are 495 NaNs out of 1,036 data points for Terra MODIS (daytime sensor; effective data points: 541); there are 369 NaNs out of 927 data points for Aqua MODIS (nighttime sensor; effective data points: 558). After removing these NaNs and combining the daytime and nighttime values (pick up the max temperature value for each day), the resulting number of the total effective data point is 744. Figure 4 displays the resulting MODIS LST time series in KST from 1999-2022 (based on 744 effective data points). An increasing trend of two-decade KST temperature time series is revealed. The temperature uptrending rate is around 0.031°C per year. In this study MODIS imagery is used to complement the other two satellite data sources despite its nature of low spatial resolution. The reason is that MODIS imagery has a shorter revisit period (sub-daily) than Landsat imagery, albeit with a coarse (1 km) spatial resolution. MODIS still provides useful information for this volcano as shown in figure 3 panel (a). In addition, all three datasets (Landsat/ASTER and MODIS) inevitably record the seasonal cycle, but MODIS shows the apparent seasonal pattern because of the denser data points. However, we have focused on the long-term trend revealed in the time series.

LST distribution on KST volcano retrieved from ASTER image dataset
The total amount of ASTER observation dataset used in this study is around 1 gigabyte. The original retrieved ASTER LST time series also contains large negative numbers which are caused by the meteorological clouds. The original retrieved ASTER LST time series contains large negative values (such as −273.14999°C; −139.24999°C; −19.94999°C, etc) that are caused by the meteorological clouds. Likewise, we have to screen out data points that are cloud contaminations marking them as NaNs. There is a total of 174 data points available for the ASTER sensor. After removing these NaNs, the resulting number of the total effective data point is 113. Figure 5 displays the resulting ASTER LST time series in KST from 1999-2022 (based on 113 effective data points). An increasing trend of two-decade KST temperature time series is revealed. The temperature uptrending rate is around 0.334°C per year.

LST distribution on KST volcano retrieved from Landsat image dataset
The total amount of Landsat Enhanced Thematic Mapper Plus (ETM+) observation dataset used in this study is around 44 gigabytes. The original retrieved Landsat LST time series also contains large negative numbers which are caused by the meteorological clouds. Likewise, we have to screen out data points that are cloud contaminations marking them as NaNs. There is a total of 190 data points available for the Landsat sensor. After removing these NaNs, the resulting number of the total effective data point is 160. Figure 6 displays the resulting Landsat LST time series in KST from 1999-2022 (based on 160 effective data points). An increasing trend of twodecade KST temperature time series is revealed. The temperature uptrending rate is around 0.038°C per year.
4. Discussion 4.1. Temporal-spatial LST distribution of KST volcano from 2000 to 2022 KST island has an area of around 2.8 km 2 , while the MODIS LST imagery datasets have a 1 km spatial resolution. As a result, MODIS sensors may be ineffective in collecting the detailed thermal activity features on KST. However, Landsat thermal infrared images (30 m pixels) are helpful to detect the thermal anomaly pattern on the land surface. Figures 7-9 shows the multitemporal LST of the KST volcano from Landsat imagery selected from 2000 to 2022. Thermally anomalous regions at the KST region of interest (ROIa; ROIb) are generally 4°C-16°C above the surroundings. Anomalous thermal distributions are indicated by the red color. The spatiotemporal analysis could also reveal whether LSTs reflect the expansion or contraction of the southern hotspot. Besides the snapshots shown in figures 7-9. Figure 10 shows the percentage of the number of LST pixels higher than 40°C at the southern KST hotspot (e.g., ROIa; total 348 pixels equivalent to 0.31 km 2 ; indicated in figure 1(a)). Variation of high LST areas (pixel counts) reflects the expansion or contraction of the southern hotspot. However, the images are generally contaminated by clouds plus gaps caused by the sensor malfunction  sources. Landsat satellite sensor holds a relatively higher differentiating ability to provide useful information on the distribution pattern of the thermal anomaly at KST, which is crucial for comprehending subsurface structures and heat flow process. In most cases, obtaining such surface anomaly patterns from ground-based point measurements is difficult due to the site's remoteness or the expense of fieldwork.

LST in the thermal anomaly region of interest ROIa and ROIb of KST volcano
To better examine the LST thermal anomaly area, we have zoomed in on the region of interest ROIa and ROIb as shown in figure 1. These two spots are determined by and extracted from the temporal-spatial LST distribution in figures 7-9, i.e., the thermal anomaly areas. Figures 11 and 12 display the Landsat LST time series in KST from 1999-2022 based on 160 effective data points. The more noticeable increasing trends of the two-decade KST Figure 9. (continued). Spatial variations in land surface temperature retrieved from Landsat imagery on KST from 2017 to 2022. Note that the images are contaminated by stripe data gaps caused by the satellite sensor malfunction. temperature time series are revealed. The temperature uptrending rates are around 0.132°C and 0.121°C per year, respectively. We shall note that the overall land and ocean temperature has risen at an average pace of 0.18°C each decade (i.e., 0.018°C per year) since 1981, according to NOAA's 2021 Annual Climate Report (NOAA 2022); nevertheless, the average rate of increase at KST in this study has been almost more than one order of magnitude that rate. In addition, according to statistics on air temperature collected over the past 20 years by Taiwan's Central Weather Bureau (CWB), The temperature ranges from 6.3°C to 34.2°C with an average of 23.5°C. This shows that the LST is primarily caused by volcanic thermal sources.

Comparison on the LST time series from MODIS, ASTER and Landsat ETM+ sensors
The magnitude of LST observation is generally reciprocal to the spatial resolution of thermal imagery. The temperature detected is the average value within the imagery pixel. Landsat and ASTER thus show higher LST temperatures than MODIS as evidenced in figures 4-6. A challenge in comparing and combining these three LST time series is the 'scale' (significant variations in magnitude and spread). Rescaling (or normalization) is needed to help us to focus on similarities of shape in each LST time series. Figure 13 demonstrates the normalization of Landsat, ASTER, and MODIS LST time series. It shows that the recorded LST temperatures from satellites are similar despite the different spatial resolutions of sensors. A long-term ground truth land surface temperature survey at KST is not available. This comparison serves as the cross-validation of the retrieved LST time series of the KST volcano.

Comparison with previous geophysical and geological surveys
Investigation on spatial and temporal distribution of low frequency volcanic earthquakes near KST shows that low frequency (LF) volcanic events were mainly located along the central graben of Southern Okinawa Trough (SOT) and to the southeast of the KST Island (Lin et al 2019). Near KST, the depths of the LF events are typically shallower than 15 km; but, as the hypocenters approach the central axis of SOT, they increase deeper. The major faults in SOT restrict the area of deep magmatic activity in addition to affecting the distribution of shallow volcanic bodies.
A survey on microearthquake activity around KST volcano shows a tight cluster of events near KST, indicating a volcano-hydrothermal system beneath KST (Konstantinou et al 2013). The majority of hypocentral depths are in the 2.5-10 km range, with the former depth equivalent to the base of the shallow sedimentary layer and the latter to the ductile lower crust. The waveforms of the three largest events demonstrate that fluid Figure 10. Percentage of the number of LST pixels whose value is higher than 40°C at KST ROIa (the southern hotspot indicated in figure 1(a)). Variation of high LST areas reflects the expansion or contraction of the southern hotspot. However, the images are generally contaminated by gaps caused by the sensor malfunction after May 31, 2003. The pixel count statistic (based on 44 data points) cannot show truthful area variations of the hotspot. The linear trend shows a decreasing (−1.188%/year) rather than an increasing tendency. movement in the upper crust causes a local stress field around KST that facilitates in the cracking along the direction of NW-SE regional extension.
A seismic study on detection of a magma reservoir beneath KST by s-wave shadows and reflections identified a km-scale molten-filled volume located beneath southeast KST (Lin et al 2018). With spatial uncertainties of 3 km, a partially-molten 19% low-velocity volume is estimated in the middle of the crust (13-23 km). The elongated orientation roughly corresponds to the strike of the Okinawa trough, suggesting that a back-arc opening could be the source of the magma reservoir.
Based on this study, the spatial LST distribution of the KST volcano indicates that LST anomaly areas are mainly located at the southeast island (as shown in figures 7-9), which is well correlated with the possible magma reservoir location estimated by the previous geophysical and geological surveys. This indicates the integration of satellite remote sensing data with geophysical and geological surveys is applicable and useful.

Satellite thermal sensor detected LSTs reflecting volcano activities underground
The land surface temperature can be contributed by the insolation, sub-surface thermal source, urban heat island (UHI) effect, local building fires, or forest fire incidents. Indeed, in most use cases of satellite thermal imagery, the temperature is often associated with changes on the Earth's surface. However, our study area is a remote non-inhabitant volcanic island, so we have ruled out other factors except the factor of the subsurface thermal source. Thermal features on volcanoes include steam vents, geysers, hot springs, lava flows, and lava domes. Satellite sensors can detect the anomalous heat released by volcanoes from space. Sometimes, before a volcanic eruption, the surface temperature at these thermal spots varies. It's possible to forecast activity by identifying these thermal anomalies. Continuous monitoring is necessary to comprehend the background / normal thermal features for detecting thermal anomalies. Infrared satellite sensors can detect volcanic thermal features, so remote sensing is a vital component for monitoring distant volcanoes with limited ground-based equipment (USGS 2022b).
Recognition of the geothermal status quo of volcanism is important for quantitatively forecasting eruption potential at volcanoes. Numerous factors make it challenging to recognize steady or non-steady state of volcanism. Previous studies (Wadge 1980, Wadge 1982) use a cumulative volume curve to assess the reliability of the steady state model for an active volcano. The cumulative volume (mass) of erupted magma must be plotted linearly with respect to time for a volcano to be in a steady state. Therefore, a consistent, well-documented eruption history and volume are needed. However, only a small number of volcanoes in the globe have received adequate study to provide this information, and the most accurate data come from observations made within the last 70 years or so. Though there is no quantitatively recorded eruption volume in KST, satellite-detected LST time series of two-decade provides a way to assess the geothermal status quo of this volcano.
In this study, the criterion for a non-steady-state at a volcano is that the plot of the LST is non-linear with respect to time as shown by the three satellite thermal sensors. In non-steady-state processes, mass and energy in the system are constantly changing with time. The conditions of heat reservoirs of KST are described by the nonsteady-state in which the heat flow that occurs while the pressures and/or rate changes with time. A system or process is in a steady state by definition if the variables (also known as state variables) that determine how it behaves are constant across time. When a flow is non-steady, its characteristics-like pressure, velocity, or density-do depend on time (Dake 1983, Chaudhry 2004).

Conclusions
The noticeable increasing trends of the two-decade KST temperature time series are revealed for three satellite thermal sensors (i.e., MODIS, ASTER, and Landsat ETM+). The calculated linear temperature uptrending rates are 0.031°C, 0.334°C, and 0.132°C per year, respectively. The upward trends shown in the two-decade LST time series of three satellite sensors signify the increasing intensity of thermal activities undergoing at the KST volcano. The spatial LST distribution of the KST volcano indicates that LST anomaly areas are mainly located on the southeast island, which is well correlated with the possible magma reservoir location from previous geophysical and geological surveys. The LST patterns imply that the activities of the subsurface thermal sources are in a non-steady-state process. In terms of methodology, this study's satellite monitoring-based approach is useful for understanding KST's eruption in the past as well as the future. The long-term temperature dataset is helpful for monitoring, hazard assessment, and mitigation of volcanic activity as well as for understanding the volcanic structure underground and geothermal activities caused by magma, gas, or fluids (e.g., water) of volcanoes. This approach is particularly useful if access to the volcano is restricted or there is a lack of funding for volcanic monitoring. Additionally, it complements and fits with the current volcanic monitoring systems, which will improve monitoring and further utilization of geothermal resources on the KST volcano.