Monitoring Thermal Activity of the Beppu Geothermal Area in Japan Using Multisource Satellite Thermal Infrared Data

: The Beppu geothermal area, one of the largest spa resorts on the northeast Kyushu Island of Japan, is fed by hydrothermal ﬂuids beneath the volcanic center of Mt. Garan and Mt. Tsurumi in the west. We explored the thermal status of the Beppu geothermal area using nighttime multisource satellite thermal infrared data (TIR) from the Advanced Spaceborne Thermal Emission and Reﬂection Radiometer (ASTER) and Landsat 8 thermal infrared scanner (TIRS) to monitor heat loss from 2009 to 2017. We also assessed heat loss from Mt. Garan fumaroles to investigate a relationship between them. The normalized differential vegetation index (NDVI) threshold method of spectral emissivity, the split-window algorithm for land surface temperature (LST), and the Stefan–Boltzmann equation for radiative heat ﬂux (RHF) were used to estimate heat loss in this study. Total heat loss increased by about a 35% trend overall from 2009 to 2015 and then declined about 33–42% in 2017 in both the Beppu geothermal area and Mt. Garan fumaroles overall. The higher thermal anomalies were found in 2015 mostly in the southeastern coastal area of the Beppu geothermal region. The highest thermal anomaly was obtained in 2011 and the lowest in 2017 within the Mt. Garan fumaroles. The areas with a higher range of RHF values were recorded in 2015 in both study areas. Finally, the results show similar patterns of heat loss and thermal anomalies in both the Beppu geothermal area and Mt. Garan fumaroles, indicating a closely connected geothermal system overall. This suggests that nighttime TIR data are effective for monitoring the thermal status of the Beppu geothermal area.


Introduction
The Beppu geothermal area is situated at the eastern edge of the Beppu-Shimabara graben that crosses from northeast to southwest on Kyushu Island in Japan ( Figure 1) [1]. Geologically, Beppu is one of the largest water-dominated geothermal systems in Japan, and it extends up from the Mt. Tsurumi-Mt. Garan volcanic center to the east coast of Beppu city, which is fed by geothermal fluid beneath the volcanic center of Mt. Tsurumi and Mt. Garan [2][3][4][5]. Beppu city is recognized historically as one of the largest spa resort areas on Kyushu Island, famous mostly for its onsen (hot spring bathing) and fumaroles [3,6]. The Beppu geothermal system is bounded by two striking faults at the north and south borders, almost along the east-west direction, and forms an alluvial fan in the center, made up by debris avalanches of andesite or dacite between the faults from the volcanoes behind [2,5]. Hydrothermal fluid flows along the two fault systems to the fan of the Beppu thermal area from the volcanic center, in southern and northern outflows known as the Beppu thermal zone

Geological Settings of the Study Area
Beppu is a tectonically controlled typical sedimentary basin at the eastern edge of the Beppu-Shimabara graben, which is under the influence of the Philippine convergence plate [23]. The Beppu-Shimabara graben extends NE-SW of the center of Kyushu Island [1]. Active normal faults are predominant in the Beppu graben along the E-W and NE-SW direction. The northern boundary of the Beppu graben is bounded by the Kannawa fault, and the southern border is bounded by the Horita and Asamagawa faults [1,6]. The Beppu geothermal field is situated within 30 km 2 on the eastern flanks of Mt. Garan and Mt. Tsurumi, which are late-Quaternary volcanic centers [8,11]. Geologically, the Beppu comprises Pliocene andesite (Kankaiji and esite), Pleistocene volcanic rocks, and sedimentary fan deposits [24,25] (Figure 2). In some parts of the northern and southern borders of the Beppu area, cretaceous granite is cropped out locally.

Geological Settings of the Study Area
Beppu is a tectonically controlled typical sedimentary basin at the eastern edge of the Beppu-Shimabara graben, which is under the influence of the Philippine convergence plate [23]. The Beppu-Shimabara graben extends NE-SW of the center of Kyushu Island [1]. Active normal faults are predominant in the Beppu graben along the E-W and NE-SW direction. The northern boundary of the Beppu graben is bounded by the Kannawa fault, and the southern border is bounded by the Horita and Asamagawa faults [1,6]. The Beppu geothermal field is situated within 30 km 2 on the eastern flanks of Mt. Garan and Mt. Tsurumi, which are late-Quaternary volcanic centers [8,11]. Geologically, the Beppu comprises Pliocene andesite (Kankaiji and esite), Pleistocene volcanic rocks, and sedimentary fan deposits [24,25] (Figure 2). In some parts of the northern and southern borders of the Beppu area, cretaceous granite is cropped out locally. Figure 2. Geological map of the study area (modified after [6]).

Satellite Data and Meteorological Data
We used daytime and nighttime visible near-infrared and thermal-infrared images of two satellite sensors, ASTER and Landsat OLI/TIRS, in this study, considering the availability of goodquality images for the first time to explore and monitor heat loss from the Beppu geothermal area and Mt. Garan fumarole of northeast Kyushu Island from 2009 to 2017. We selected three pairs (day and night) of ASTER and two pairs (day and night) of Landsat 8 OLI/TIRS images considering cloudfree, good-quality, and same-season images of the study area ( Table 1). The images were best quality (ranked 9 out of 9) (https://landsat.usgs.gov/landsat-8-l8-data-users-handbook-section-4) and 100% cloud-free, and were acquired from the United States Geological Survey (USGS) archives that were both radiometrically and geometrically corrected. Nighttime thermal infrared images were used in this study to avoid or reduce the solar effect to calculate geothermal heat loss from the volcanic thermal grounds. The daytime visible near-infrared images were used to retrieve the emissivity for

Satellite Data and Meteorological Data
We used daytime and nighttime visible near-infrared and thermal-infrared images of two satellite sensors, ASTER and Landsat OLI/TIRS, in this study, considering the availability of good-quality images for the first time to explore and monitor heat loss from the Beppu geothermal area and Mt. Garan fumarole of northeast Kyushu Island from 2009 to 2017. We selected three pairs (day and night) of ASTER and two pairs (day and night) of Landsat 8 OLI/TIRS images considering cloud-free, good-quality, and same-season images of the study area (Table 1). The images were best quality (ranked 9 out of 9) (https://landsat.usgs.gov/landsat-8-l8-data-users-handbook-section-4) and 100% cloud-free, and were acquired from the United States Geological Survey (USGS) archives that were both radiometrically and geometrically corrected. Nighttime thermal infrared images were used in this study to avoid or reduce the solar effect to calculate geothermal heat loss from the volcanic thermal grounds. The daytime visible near-infrared images were used to retrieve the emissivity for each land cover of the study area, as emissivity is one of the important factors for heat loss retrieval from thermal ground. The ASTER satellite has 14 multispectral bands: three visible near-infrared bands of 15 m, six shortwave infrared bands of 30 m, and five bands of 90 m spatial resolution. The latest Landsat 8 satellite carries OLI and TIRS sensors, ensuring continuity of image collection, which have nine OLI bands of 30 m image resolution and three TIRS bands of 100 m resolution. We also acquired ASTER standard products (Emissivity, AST_05; Surface kinetic temperature, AST_08) from the land processes distributed active archive center (LP DAAC) of USGS for validating the heat loss monitoring from the study areas. The meteorological data, such as relative humidity and ambient temperature, were collected from the Japanese Automated Meteorological Data Acquisition System (AMeDAS) Oita station, and these data were later adjusted with altitude variation between the station (4.6 m above mean sea level) and the Mt. Garandake fumaroles (1045 m above mean sea level) and the Beppu thermal area using the international standard atmospheric (ISA) method [6,26,27] (Figure 1). Path/Row is the image acquisition system of the world by satellite as described for Landsat at: https://landsat.gsfc. nasa.gov/the-worldwide-reference-system/.

Emissivity Retrieval from ASTER Data
Both sensors of satellite images were analyzed for heat flow measurements by using various steps after collecting them from the USGS archives free of cost ( Figure 3). ASTER daytime images were used to estimate spectral emissivity using the normalized differential vegetation index (NVDI) threshold method for the study areas first, so the necessary atmospheric correction and conversion of digital number (DN) value to reflectance was done for the bands of red and near-infrared, and these reflectance values were later used to calculate NDVI values of the study area for the respective year's images. The NDVI value was used for estimating emissivity of the study area from 2009 to 2017. The spectral emissivity of each land cover of the study areas was estimated using the NDVI threshold method according to the following equation for thermal bands 13 and 14 of ASTER data [28]: where A v = (1 − A s ), and A v and A s are the vegetation and soil coverage of any pixel; A v = (NDVI − NDVI s )/(NDVI v − NDVI s ), and NDVI v and NDVI s are NDVI of vegetation and soil, respectively; NDVI = (b 3 − b 2 )/(b 3 + b 2 ), and b 3 and b 2 are the infrared band and red band of ASTER data, respectively. In the normal case, the NDVI value is above 0.75 for full healthy vegetation coverage of a pixel and below 0.2 for full bare ground or soil. R v = 0.92762 + 0.07033A v and R s = 0.99782 + 0.08362A v . ε s13 and ε s14 are the spectral emissivity of soil for bands 13 and 14, respectively (ε s13 = 0.9676 and ε s14 = 0.9779), and ε v13 and ε v14 are the spectral emissivity of vegetation for the bands 13 and 14, respectively (ε s13 = 0.9867 and ε s14 = 0.9899).  Figure 3. Satellite image processing flowchart of the study. Both day-and nighttime images from ASTER and Landsat 8 sensors were used to estimate spectral emissivity and land surface temperature using the various methods listed above. NIR, near-infrared; TIR, thermal infrared; NDVI, normalized differential vegetation index; RHF, radiative heat flux; HDR, heat discharge rate.

LST Retrieval from ASTER Data
We used the nighttime ASTER thermal infrared data to estimate true land surface temperature of the study areas using the split-window algorithm, taking into consideration thermal bands 13 and 14 of ASTER images. First, the thermal band data were used to convert into radiance values from the DN of the selected nighttime satellite images using the metadata of those images with the necessary equation [28]. The radiance value was used to estimate the brightness temperature of the study area from the sensor's information without considering the atmospheric effects. Then, the brightness temperature was used to estimate true skin temperature of the ground using the following splitwindow algorithm of [29] with TIR bands 13 and 14 of ASTER images: , εi and τi are the emissivity atmospheric transmissivity of ASTER band i. τ13 = 0.979160 − 0.062918w and τ14 = 0.968144 − 0.098942w (in case of water vapor content (w) 0.4-2.0). The total water vapor content was calculated using the relative air humidity, air density, and air saturation mix ratio from local meteorological data [28].

Emissivity Retrieval from Landsat 8 OLI Data
For the Landsat 8 OLI-TIRS images, we used daytime OLI images to estimate the spectral emissivity of each land cover of the study area and nighttime TIRS images to estimate the land surface temperature without reducing the solar effect. First the digital values of the visible near-infrared . Satellite image processing flowchart of the study. Both day-and nighttime images from ASTER and Landsat 8 sensors were used to estimate spectral emissivity and land surface temperature using the various methods listed above. NIR, near-infrared; TIR, thermal infrared; NDVI, normalized differential vegetation index; RHF, radiative heat flux; HDR, heat discharge rate.

LST Retrieval from ASTER Data
We used the nighttime ASTER thermal infrared data to estimate true land surface temperature of the study areas using the split-window algorithm, taking into consideration thermal bands 13 and 14 of ASTER images. First, the thermal band data were used to convert into radiance values from the DN of the selected nighttime satellite images using the metadata of those images with the necessary equation [28]. The radiance value was used to estimate the brightness temperature of the study area from the sensor's information without considering the atmospheric effects. Then, the brightness temperature was used to estimate true skin temperature of the ground using the following split-window algorithm of [29] with TIR bands 13 and 14 of ASTER images: where T 13  . The total water vapor content was calculated using the relative air humidity, air density, and air saturation mix ratio from local meteorological data [28].

Emissivity Retrieval from Landsat 8 OLI Data
For the Landsat 8 OLI-TIRS images, we used daytime OLI images to estimate the spectral emissivity of each land cover of the study area and nighttime TIRS images to estimate the land surface temperature without reducing the solar effect. First the digital values of the visible near-infrared bands of daytime Landsat 8 OLI images were converted into reflectances with necessary information from metadata, then the spectral emissivity of the land covers of the study area was estimated using the NDVI threshold method for thermal bands 10 and 11 of Landsat TIRS images as follows [30]: where NDVI = ( nir − r )/( nir + r ), here, nir and r are the reflectance of near-infrared and red band; P v is the vegetation fraction (NDVI − NDVI min /NDVI max − NDVI min ) 2 , NDVI max = 0.5, ; and F' is the geometric factor ranging between 0 and 1 (typically 0.55) [30]. For bare soil (NDVI < 0.2 and P v = 0), the emissivity of bands 10 and 11 was ε 10 = 0.973 − 0.047δ 4 and ε 11 = 0.984 − 0.026δ 4 , and δ i is the reflectance of the Landsat 8 OLI bands [30].

LST Retrieval from Landsat 8 TIRS Data
We used nighttime Landsat 8 thermal infrared data to estimate the land surface temperature, which was later used for geothermal heat flux measurement from the Garandake-Beppu thermal area with or without reducing solar effect. First thermal infrared bands 10 and 11 of Landsat 8 were converted into radiance using the equation L λ = M L × Q cal + A L ; here, L λ is the spectral radiance (W/m 2 ·sr·µm), M L is the multiplicative scaling factor, Q cal is the pixel value, and A L is the radiance additive factor (taken from the metadata of the image). Then the radiance value was converted into effective or brightness temperature using the equation T i = K 2 /In(K 1 /L λ + 1); here, T i is the brightness temperature, K 1 and K 2 are thermal efficient coefficients (taken from the metadata), and L λ is the radiance value. We used local meteorological data from the Oita AMeDAS station, such as ambient temperature and relative humidity, to calculate atmospheric transmissivity. Because of many difficulties in obtaining in situ atmospheric transmissivity values when the satellite passes, we used simulated equations to get atmospheric transmissivity using the water vapor content in the following equation [30]: where τ i is the atmospheric transmissivity of the respected bands of the Landsat thermal infrared data, w is total water vapor content (g/m 2 ), H is relative air humidity (%), E is the saturation mix ratio (g/Kg) of water vapor, A is air density (g/m 3 ), and R w (0) = 0.6834 and 0.6356 for mid-latitude summer and winter (Table 2). Table 2. Relationship between air temperature (T), saturation mix ratio (E), and air density (A) [30]. We used a split-window algorithm to estimate skin LST of the study areas using Landsat 8 TIRS bands 10 and 11, which was developed based on Qin et al., 2001, as follows [30]: where T s is the LST; , ε is emissivity and τ is atmospheric transmissivity; and Li is the linear fitting coefficient taken from [30].

RHF and Heat Loss Retrieval
Actually, the total heat loss or heat discharge rate (HDR) is the summation of conductive, convective, and radiative heat loss from any volcano or geothermal field without solar effect. As satellite thermal infrared data can only be used to study radiative heat loss, radiative heat flux (RHF) was calculated using the Stefan-Boltzmann equation within the Garandake-Beppu thermal area from 2009 to 2017 [14,31,32]. Strictly speaking, RHF can be defined as the radiative portion of surface heat loss without solar effect from any volcano or geothermal field. The Stefan-Boltzmann equation defines that the thermal energy radiated by a true blackbody per second per unit area is proportional to the fourth power of absolute temperature. According to the Stefan-Boltzmann law, if the hot object or geothermal ground is other than an ideal radiator or blackbody and the thermal surface is radiating to its cooler surrounding atmosphere at temperature T a , then the net RHF from the thermal area would be as follows [31,32]: where Q is the radiative heat flux, σ is the Stefan-Boltzmann constant (5.6703 × 10 −8 watt/m 2 K 4 ), A is the area of radiation, T s is the surface or hot body temperature, and T a is the cooler surrounding temperature. We calculated total radiative heat loss (RHL) by summing the pixels' positive RHF values for each year's images of the study area at this stage. We used a relationship coefficient (about 15%) between radiative heat flux and heat discharge rate to calculate total heat loss and monitor the trend of heat discharge rate in the study area [16].

Results and Discussion
Nighttime multispectral satellite thermal infrared data are used for the first time in this study to investigate and monitor the thermal status of the Beppu geothermal area and the adjacent Mt. Garan (Garandake) fumaroles of the northeast part of Kyushu Island from 2009 to 2017. The intention of nighttime thermal infrared data is to reduce or remove the solar effect to monitor actual geothermal heat loss in the studied area. Daytime images from the same month as nighttime images were used to estimate the spectral emissivity for land cover of the study areas. As the Beppu geothermal system is proximal to the Mt. Garan and Mt. Tsurumi volcano in the west, we also assessed the thermal heat loss of Mt. Garan fumaroles to understand the relationship between them. As the elevation variation between the Beppu geothermal area (less than 300 m) and Mt. Garan fumarole (above 1000 m) is large, we split the study area in two. The results are discussed separately to better understand the relationship between them.

Mt. Garan Fumaroles
Five daytime visible near-infrared (VNIR) images of the Garandake fumaroles were acquired from USGS archives, three ASTER and two Landsat 8 OLI, considered to be good quality and cloud-free ( Figure 4). Daytime images were used to calculate the spectral emissivity of each land cover of the study area using the NDVI threshold method. NDVI indicates the ratio of near-infrared and red bands of ASTER multispectral satellite images. NDVI values range from −1 to 1. Higher NDVI values (above 0.5) indicate healthy vegetation, lower values (below 0.2) indicate bare land or water bodies and mixed land cover (NDVI = 0.2-0.5) [32]. The Garandake fumaroles showed lower or negative NDVI values, indicating bare ground around the fumaroles, and higher positive values in vegetated areas ( Figure 5). The NDVI threshold method were used to retrieve spectral emissivity of each land cover of the study area. Average spectral emissivity ranges from 0.988 to 0.94 for the land covers of the Mt. Garan fumaroles area ( Figure 6).
Land surface temperature is the main component of heat loss measurement from volcanoes or volcanically thermal ground. Higher spatial resolution of the thermal infrared band is the best option for detailed or accurate results of thermal anomalies or heat loss component. ASTER or Landsat 8 thermal infrared data consist of moderate-resolution images (resampled 30 m of original 90-100 m TIR data) and are used nowadays effectively for the thermal mapping of volcanoes. We estimated the thermal features, i.e., land surface temperatures, of the study areas with the ASTER and Landsat 8 TIRS data using a split-window algorithm with two thermal bands of ASTER (bands 13 and 14) and Landsat 8 TIRS (bands 10 and 11) images. We obtained a higher LST anomaly without ambient temperature within the Garandake fumaroles area than the surrounding area in all five sets of thermal images from 2009 to 2017. The highest LST without ambient temperature in the Garandake fumaroles was about 11 • C in 2011 and the lowest was about 7 • C in 2017 ( Figure 7; Table 3). The thermal anomaly area was also the highest in 2011 and the lowest in 2017. TIR data) and are used nowadays effectively for the thermal mapping of volcanoes. We estimated the thermal features, i.e., land surface temperatures, of the study areas with the ASTER and Landsat 8 TIRS data using a split-window algorithm with two thermal bands of ASTER (bands 13 and 14) and Landsat 8 TIRS (bands 10 and 11) images. We obtained a higher LST anomaly without ambient temperature within the Garandake fumaroles area than the surrounding area in all five sets of thermal images from 2009 to 2017. The highest LST without ambient temperature in the Garandake fumaroles was about 11 °C in 2011 and the lowest was about 7 °C in 2017 ( Figure 7; Table 3). The thermal anomaly area was also the highest in 2011 and the lowest in 2017.         Radiative heat flux is one of the components of heat loss from volcanoes that can be measured using satellite thermal infrared sensors. We applied the Stefan-Boltzmann law of thermal heat loss to estimate RHF using skin land surface temperature of the Garandake fumaroles from 2009 to 2017 derived from the ASTER and Landsat 8 thermal infrared data. The highest pixel RHF value was within the range of 47 to 53 W/m 2 from 2009 to 2015, and 38 W/m 2 in 2017 ( Figure 8; Table 3 (Table  3). Otherwise, we obtained RHL of about 16, 12, and 12 MW in 2011, 2012, and 2017, respectively ( Figure 9). Then the total heat loss or heat discharge rate (HDR) of the Garandake fumarolic area was estimated after multiplying total RHL by a relationship coefficient between RHF and HDR. We  Radiative heat flux is one of the components of heat loss from volcanoes that can be measured using satellite thermal infrared sensors. We applied the Stefan-Boltzmann law of thermal heat loss to estimate RHF using skin land surface temperature of the Garandake fumaroles from 2009 to 2017 derived from the ASTER and Landsat 8 thermal infrared data. The highest pixel RHF value was within the range of 47 to 53 W/m 2 from 2009 to 2015, and 38 W/m 2 in 2017 ( Figure 8; Table 3 (Table 3). Otherwise, we obtained RHL of about 16, 12, and 12 MW in 2011, 2012, and 2017, respectively ( Figure 9). Then the total heat loss or heat discharge rate (HDR) of the Garandake fumarolic area was estimated after multiplying total RHL by a relationship coefficient between RHF and HDR. We obtained the highest heat loss of about 112 MW in 2015 and the lowest of about 73 MW in 2009 (Table 3; Figure 9). There were two consecutive peaks of heat loss from the Mt. Garan fumaroles, one in 2011 and the other in 2015 (Figure 9). In the higher ranges of RHF anomaly, the maximum was in 2015 and the lowest in 2017; the mid-ranges of RHF anomaly were the highest in 2015 and lowest in 2009, and the lower range of RHF anomaly were at the maximum in 2017 and minimum in 2015 ( Figure 10).

Beppu Geothermal Area
The Beppu thermal area covers almost all of Beppu city. A false color composite map of the Beppu thermal area shows vegetation as red, urban areas as bluish green, and water bodies as blue for ASTER and Landsat OLI images from 2009 to 2017 ( Figure 11). Spectral emissivity is an important component to be considered for heat flux estimation in any volcanic or thermal area. We estimated the spectral emissivity of the Beppu thermal area using the NDVI threshold method for respective thermal bands of ASTER and Landsat 8. NDVI is derived using the red and near-infrared bands of daytime images of ASTER and Landsat OLI images of the Beppu thermal area from 2009 to 2017 ( Figure 12). NDVI values range from −1 to +1. A higher positive NDVI indicates healthy vegetation, and negative or close to zero indicates bare or urban area or mixed land. We obtained an average spectral emissivity of the Beppu geothermal area from 0.93 to 0.99 ( Figure 13). The higher value indicates healthy vegetation and the lower value indicates bare or urban or mixed land of the Beppu geothermal area.
Land surface temperature was estimated after applying the split-window algorithm with the nighttime thermal infrared data of ASTER and Landsat 8 TIRS images of the Beppu geothermal area  (Table  4). Higher RHF anomaly was observed in the southeastern part of the Beppu thermal area ( Figure  15). Total radiative heat loss was estimated after summation of all positive RHF pixels.  Table 4). The highest total heat loss or heat discharge rate was in 2015, about 1279 MW, and the lowest was about 739 MW in 2017 ( Figure  16). Dividing the RHF anomaly areas into three parts, high, moderate, and low, we found that the high heat loss anomaly area converged to maximum in 2015, the moderate anomaly area in 2017, and the lower RHF anomaly area in 2011 ( Figure 17).

Beppu Geothermal Area
The Beppu thermal area covers almost all of Beppu city. A false color composite map of the Beppu thermal area shows vegetation as red, urban areas as bluish green, and water bodies as blue for ASTER and Landsat OLI images from 2009 to 2017 ( Figure 11). Spectral emissivity is an important component to be considered for heat flux estimation in any volcanic or thermal area. We estimated the spectral emissivity of the Beppu thermal area using the NDVI threshold method for respective thermal bands of ASTER and Landsat 8. NDVI is derived using the red and near-infrared bands of daytime images of ASTER and Landsat OLI images of the Beppu thermal area from 2009 to 2017 ( Figure 12). NDVI values range from −1 to +1. A higher positive NDVI indicates healthy vegetation, and negative or close to zero indicates bare or urban area or mixed land. We obtained an average spectral emissivity of the Beppu geothermal area from 0.93 to 0.99 ( Figure 13). The higher value indicates healthy vegetation and the lower value indicates bare or urban or mixed land of the Beppu geothermal area.
Land surface temperature was estimated after applying the split-window algorithm with the nighttime thermal infrared data of ASTER and  Table 4). The highest total heat loss or heat discharge rate was in 2015, about 1279 MW, and the lowest was about 739 MW in 2017 ( Figure 16). Dividing the RHF anomaly areas into three parts, high, moderate, and low, we found that the high heat loss anomaly area converged to maximum in 2015, the moderate anomaly area in 2017, and the lower RHF anomaly area in 2011 ( Figure 17).

Validation of the Results
ASTER standard products were used to validate in this study, just for the results of LST and RHL, retrieved using ASTER images for both study areas. We used nighttime high level products of emissivity (AST_05) and surface kinetic temperature (AST_08) to retrieve LST and RHF for the year of 2009, 2011 and 2013 to validate the results of observed heat loss from both Mt. Garan fumaroles and Beppu geothermal area in this study ( Figure 18). Although the anomalies of LST and RHF are not closely matched between the results of observed and standard products of ASTER (may be due to the resolution difference 15 and 90 m respectively for observed and standard), but the total RHL showed similar pattern and close to the observed total RHL of the study areas (Figures 18-20) (Table  5). We analyzed the Pearson correlation coefficient between total RHL based on ASTER standard products and observed total RHL, and obtained strong correlation (0.9822 for Mt. Garan and 0.9881 for the Beppu geothermal area) between them for the study areas ( Figure 21; Table 5). The variation between total RHL based on ASTER standard products and observed total RHL may be related with the atmospheric profile and altitude differences. We used the local Oita meteorological station's data (AMeDAS, used to acquire atmospheric data in every 10 min) with the same time of the image acquisition in our study, but the meteorological station is located about 15-20 km southeast of the study area. We were also adjusted the altitude variation of atmospheric ambient temperature for the whole Beppu area using a single altitude (284 m as average) for the whole Beppu region. Hence, the total RHLs are varied more between observed and ASTER standard product's base RHL in the Beppu geothermal area (large area with varied topography) than Mt. Garan fumaroles (small area and more or less uniform topography) ( Figure 20).

Validation of the Results
ASTER standard products were used to validate in this study, just for the results of LST and RHL, retrieved using ASTER images for both study areas. We used nighttime high level products of emissivity (AST_05) and surface kinetic temperature (AST_08) to retrieve LST and RHF for the year of 2009, 2011 and 2013 to validate the results of observed heat loss from both Mt. Garan fumaroles and Beppu geothermal area in this study ( Figure 18). Although the anomalies of LST and RHF are not closely matched between the results of observed and standard products of ASTER (may be due to the resolution difference 15 and 90 m respectively for observed and standard), but the total RHL showed similar pattern and close to the observed total RHL of the study areas (Figures 18-20) (Table 5). We analyzed the Pearson correlation coefficient between total RHL based on ASTER standard products and observed total RHL, and obtained strong correlation (0.9822 for Mt. Garan and 0.9881 for the Beppu geothermal area) between them for the study areas ( Figure 21; Table 5). The variation between total RHL based on ASTER standard products and observed total RHL may be related with the atmospheric profile and altitude differences. We used the local Oita meteorological station's data (AMeDAS, used to acquire atmospheric data in every 10 min) with the same time of the image acquisition in our study, but the meteorological station is located about 15-20 km southeast of the study area. We were also adjusted the altitude variation of atmospheric ambient temperature for the whole Beppu area using a single altitude (284 m as average) for the whole Beppu region. Hence, the total RHLs are varied more between observed and ASTER standard product's base RHL in the Beppu geothermal area (large area with varied topography) than Mt. Garan fumaroles (small area and more or less uniform topography) (Figure 20).

Discussion
The results of this research could be more effective to monitor the thermal activity of the Mt. Garan fumaroles and Beppu geothermal area if it was possible to obtain same sensor of satellite nighttime images of each year during the study period. We did not consider the daytime thermal infrared data due to the solar insolation as well as reduce or remove the solar effect to monitor the geothermal activity of the study area. To get simultaneous images of both sensors of any specific year could be another finding for sensor variations of heat loss from the study area but it was not possible to acquire like these images for a year. We were just able to get good quality and cloud-free ASTER images only for the year of 2009, 2011 and 2013, and Landsat images for the year of 2015, 2017. Although the resolution of thermal infrared data was 90 m for ASTER and 100 m for Landsat 8 TIRS, but we were able to acquire with same resolution of the thermal infrared data as 30 m (resampled) from the USGS archive in case of both sensors.
It is well known that LST is the important component for geothermal heat loss calculation that can be retrieve using the thermal infrared data of various sensors such as Landsat, ASTER among others. There are several types of algorithm nowadays to retrieve LST, among them, the temperature and emissivity separation (TES) and split-window (SW) algorithm are widely used to retrieve LST using ASTER thermal infrared image. It is recognized that the TES algorithm is designed for the ASTER image to calculate LST and emissivity using a maximum-minimum difference (MMD) module, which is an empirical relationship between the minimum emissivity and the difference between maximum and minimum emissivity, to make the number of observe equations equal to that of the unknown variables [33]. This method can estimate LST and emissivity only with atmospheric-corrected data of TIR without any prior knowledge of emissivity [33]. The accuracy of atmospheric correction is an important factor for the efficiency of the TES method, and sometimes, the results of this algorithm are imprecise for low emissivity land cover such as dense vegetation, water and snow [33,34]. On the other hand, the SW algorithm can estimate LST by removing atmospheric effects from the linear or nonlinear combination of the brightness temperature of the nearby bands (i.e., like 11 and 12 µm in wavelength) and therefore, reduces the requirements of detail atmospheric data [33]. The SW method is used widely nowadays, because of its simplicity, efficiency and insensitive to atmospheric uncertainty from a variety of thermal infrared sensors such as Landsat, ASTER among others [13,14,30]. The only prior requisite of SW method is the emissivity of land cover in advance of the study area. This is the motivation to use the SW method for LST estimation in this study using ASTER thermal infrared data. As the land surface emissivity is one of the key component to retrieve LST from remote sensing data. There are various efficient methods to retrieve land surface emissivity from normalized differential vegetation index (NDVI) values [35,36]. The NDVI threshold based emissivity method is very famous and was used in this study to retrieve emissivity for soil, vegetation and water bodies as 15 m in resolution for ASTER and 30 m resolution for Landsat 8 images with efficiently [13,14,28,36], contrasting the TES method used to retrieve emissivity of 90 m in resolution. In case of Landsat 8 OLI/TIRS sensor, it is documented that the TIRS suffered from stray light problem since launch of this satellite [37]. A developed algorithm was applied to correct the stray light problem and supplied the corrected TIRS data from February 2017, but it was not recommended to apply the TIRS band 11 for split-window algorithm [37]. With this in mind, according to our intention to use same method for both ASTER and Landsat 8 TIRS sensor, we applied the SW algorithm using two thermal bands of both sensors in this study from 2009 to 2017.
Geothermal heat loss monitored using nighttime ASTER and Landsat 8 TIRS data from 2009 to 2017 shows that the highest total heat loss occurred in 2015 and the lowest in 2017, although the highest maximum pixel RHF was recorded in 2015 for both the Beppu geothermal area and Mt. Garan fumaroles, and the lowest in 2011 for the Beppu geothermal area and in 2017 for the Mt. Garan fumaroles (Tables 3 and 4; Figures 9 and 16). Thus, thermal activity increased from 2009 to 2015 in both the Beppu thermal area and Mt. Garan fumaroles and declined afterward up to 2017. We obtained a moderate relationship coefficient (0.59) between the observed total RHL of Mt. Garan fumaroles and Beppu geothermal area using the Pearson correlation coefficient analysis free software of [38] ( Figure 22). Spatial analysis of RHF showed that maximum values of higher and middle ranges of RHF were observed in 2015, and lower ranges in 2017 in the Mt. Garan fumaroles ( Figure 10). On the other hand, the high heat loss anomaly area converged to maximum in 2015, moderate anomaly area in 2017, and lower RHF anomaly area in 2011 in the Beppu geothermal area (Figure 17). Nighttime image analysis showed that the highest thermal anomaly was associated with the Beppu thermal area, mostly along southeast coastal region, although part of the northern region around the Chinoike Jigoku showed high value, with maximum recorded LST without ambient temperature about 7 • C in 2015 and lowest about 5 • C in 2017 ( Figure 14). On the other hand, thermal anomaly showed the highest maximum LST of about 11 • C in 2011 without ambient temperature and lowest of about 7 • C in 2017 within the fumarolic area of Mt. Garan (Figure 7). We obtained spectral emissivity in the range of 0.93 to 0.988 in the Beppu geothermal area, and from 0.94 to 0.988 in the Mt. Garan fumaroles area, where the vegetation showed higher values and bare ground or water bodies showed lower values generally (Tables 3 and 4). Because of the sensors and methods variation between ASTER and Landsat 8, the results of NDVI and emissivity may vary little bit for the land covers of the study area. Understanding the relationship between RHF and LST against NDVI (or land cover), we observed that the NDVI values less than 0.5 (mixed or bare land) were associated with highest LST and RHF of the Beppu geothermal area and Mt. Garan fumaroles. The study indicates that the Beppu geothermal area is closely connected with the Mt. Garan fumaroles area on the basis of the trend of total heat loss and thermal anomalies from 2009 to 2017.
There are some limitations of this study. A significant limitation of this research is the lack of ground validation due to the inaccessibility of the fumaroles, but we validated the results of ASTER image retrieved total RHL for both of our study regions using the ASTER standard products (AST_05, emissivity and AST_08, surface kinetic temperature) retrieve total RHL in this study. We were also obtained similar results of RHL (11)(12)(13)(14)(15)(16)(17) in a recent study of airborne thermal infrared data based heat loss (20 MW) from the Mt. Garan fumaroles in 1987 [12]. Another limitation was that we did not consider urban heat loss in this study for the city in the Beppu geothermal area, although we used nighttime satellite thermal infrared data to reduce the effect of urbanization and solar insolation. With all these limitations, as a first study for thermal activity monitoring using satellite images, we recommend that the applied method would be worthy for long term monitoring of the thermal activity using continuous satellite nighttime thermal infrared data for the sustainability of the resources and utilization within the Mt. Garan fumaroles and the Beppu geothermal area.  (Tables 3 and 4). Because of the sensors and methods variation between ASTER and Landsat 8, the results of NDVI and emissivity may vary little bit for the land covers of the study area. Understanding the relationship between RHF and LST against NDVI (or land cover), we observed that the NDVI values less than 0.5 (mixed or bare land) were associated with highest LST and RHF of the Beppu geothermal area and Mt. Garan fumaroles. The study indicates that the Beppu geothermal area is closely connected with the Mt. Garan fumaroles area on the basis of the trend of total heat loss and thermal anomalies from 2009 to 2017. There are some limitations of this study. A significant limitation of this research is the lack of ground validation due to the inaccessibility of the fumaroles, but we validated the results of ASTER image retrieved total RHL for both of our study regions using the ASTER standard products (AST_05, emissivity and AST_08, surface kinetic temperature) retrieve total RHL in this study. We were also obtained similar results of RHL (11)(12)(13)(14)(15)(16)(17) in a recent study of airborne thermal infrared data based heat loss (20 MW) from the Mt. Garan fumaroles in 1987 [12]. Another limitation was that we did not consider urban heat loss in this study for the city in the Beppu geothermal area, although we used nighttime satellite thermal infrared data to reduce the effect of urbanization and solar insolation. With all these limitations, as a first study for thermal activity monitoring using satellite images, we recommend that the applied method would be worthy for long term monitoring of the thermal activity using continuous satellite nighttime thermal infrared data for the sustainability of the resources and utilization within the Mt. Garan fumaroles and the Beppu geothermal area.