Heat transport process associated with the 2021 eruption of Aso volcano revealed by thermal and gas monitoring

The thermal activity of a magmatic–hydrothermal system commonly changes at various stages of volcanic activity. Few studies have provided an entire picture of the thermal activity of such a system over an eruptive cycle, which is essential for understanding the subsurface heat transport process that culminates in an eruption. This study quantitatively evaluated a sequence of thermal activity associated with two phreatic eruptions in 2021 at Aso volcano. We estimated plume-laden heat discharge rates and corresponding H 2 O flux during 2020–2022 by using two simple methods. We then validated the estimated H 2 O flux by comparison with volcanic gas monitoring results. Our results showed that the heat discharge rate varied substantially throughout the eruptive cycle. During the pre-eruptive quiescent period (June 2020–May 2021), anomalously large heat discharge (300–800 MW) were observed that were likely due to enhanced magma convection degassing. During the run-up period (June–October 2021), there was no evident change in heat discharge (300–500 MW), but this was accompanied by simultaneous pressurization and heating of an underlying hydrothermal system. These signals imply progress of partial sealing of the hydrothermal system. In the co-eruptive period, the subsequent heat supply from a magmatic region resulted in additional pressurization, which led to the first eruption (October 14, 2021). The heat discharge rates peaked (2000–4000 MW) the day before the second eruption (October 19, 2021), which was accompanied by sustained pressurization of the magma chamber that eventually resulted in a more explosive eruption. In the post-eruptive period, enhanced heat discharge (~ 1000 MW) continued for four months, and finally returned to the background level of the quiescent period (< 300 MW) in early March 2022. Despite using simple models, we quantitatively tracked transient thermal activity and revealed the underlying heat transport processes throughout the Aso 2021 eruptive activity.


Introduction
Thermal activity of active volcanoes persists regardless of whether they are in an eruptive or quiescent period (Kagiyama 1981;Rowe et al. 1992;Jousset et al. 2000).Heat discharges during eruptive and quiescent periods are on the same order of magnitude (~ 100 MW/100 km in Japan Arc; Kagiyama 1986), highlighting the importance of quantifying heat discharge equally for both periods.Enhanced thermal activity has been observed during transitional periods from quiescence to magmatic eruption (Dehn et al. 2002;Hernandez et al. 2007;Wessels et al. 2010).This is also the case for phreatic eruptions, for which precursory thermal manifestations are relatively unclear, with small magnitudes and various leading times (Barberi et al. 1992).
Recent studies have successfully captured slight thermal precursors of phreatic eruptions.However, in many cases only qualitative records of thermal activity sequences are provided, such as temperature increases in crater lakes, springs, and fumaroles (Christenson et al. 2010;Rouwet et al. 2021) and the geothermal anomaly area expansions (Tajima et al. 2020).Several studies have attempted to quantify the heat discharge at different stages, such as the net change before and after an eruption (Ehara et al. 2005;Mannen et al. 2018), co-eruptive heat discharge (Terada et al. 2021), post-eruptive temporal change (Nakaboh et al. 2003;Narita et al. 2019), and rapid increases due to unrest events (Chiodini et al. 2001;Inguaggiato et al. 2012;Behr et al. 2023).However, these studies obtained only snapshots of the thermal activity over an eruptive cycle (i.e., a period spanning before, during, and after an eruption).Only a few studies have quantified the escalation of pre-eruptive thermal activity to provide a complete picture of the thermal activity over an eruptive cycle (Aguilera et al. 2021;Thompson et al. 2022).
In magmatic-hydrothermal systems, which are responsible for phreatic eruptions, thermal energy is discharged in various forms, such as fumarole fields (Inguaggiato et al. 2023), geothermal anomalies (Harris et al. 2009), gas plumes (Matsushima et al. 2003), and hot water ponds and crater lakes (Hurst et al. 1991;Tajima et al. 2020).In recent years, satellite-based thermal infrared (TIR) imagery has become a powerful tool for monitoring such thermal phenomena in hydrothermal systems (Gaete et al. 2020;Coppola et al. 2022).However, this tool does not provide information on the gas-laden heat discharge.In persistently degassing volcanoes, a large portion of heat is discharged as gas plumes (e.g., Kagiyama 1981).Thus, accurate estimation of the plume heat discharge is necessary to understand the thermal process that culminates in eruptions.
At the Nakadake crater of Aso volcano, south of Japan, plume discharge continued from 2020 to 2022, interposed by phreatic eruptions in 2021.Previous studies on the Nakadake crater have investigated the thermal activity of plume (Kagiyama 1981;Fukui 1995) and a crater lake (Saito et al. 2008;Terada et al. 2008) during quiescent periods and fumaroles over the southern wall of the crater (Yokoo and Ishii 2021).However, the pre-eruptive thermal precursors have not yet been thoroughly investigated at Aso volcano.Cigolini et al. (2018) analyzed satellite TIR images over the crater during the 2014-2016 eruption, but they did not provide absolute values of the plume heat discharge.Nashimoto and Yokoo (2023) analyzed aerial TIR imagery in 2020-2022, but they only targeted the fumarolic fields.Considering active plume discharge, these plumes were probably the primary form of heat discharge at the crater during 2020-2022.Thus, accurately estimating their contribution is required to obtain the total heat discharge by the magmatic-hydrothermal system of Aso volcano.
In this study, we tracked the temporal evolution of thermal activity in the Nakadake crater in 2020-2022 to reveal the underlying mechanism of heat transport responsible for two phreatic eruptions in 2021.Our observations targeted the plume emissions from vents in the crater floor.We estimate the plume heat discharge rate (HDR) based on the buoyant plume-rise and TIR image-based models.The use of multiple methods to monitor long-term thermal activity (i.e. on the scale of years) has not yet been fully addressed.We validated our estimates by comparing them with measurements of the SO 2 gas flux and gas composition.Finally, we discuss the association between surface thermal activity and the 2021 eruptions based on the temporal evolution of the HDR and other geophysical observations.

Recent activity of Aso volcano
Aso caldera formed from four caldera eruptions during 266-89 ka, with post-caldera activity concentrated at the center of the caldera (Ono and Watanabe 1985;Ono et al. 1995;Shinmura et al. 2022) (Fig. 1a).The Nakadake crater has been the most active crater and center of recent thermal activity and eruptions.During quiescent periods, the crater hosts a hot acidic lake called 'Yudamari' , of which water level changes in response to volcanic activity (Ikebe and Watanabe 1990;Terada et al. 2008).Recent cycles of magmatic activity began in 2014; Strombolian eruptions and continuous magmatic ash emissions, locally called ' Ash eruption' at Aso volcano (Ono et al. 1995), continued from November 2014 to May 2015 (Miyabuchi and Hara 2019;Ishii and Yokoo 2021).The activity then shifted to phreatomagmatic eruptions in September 2015 and October 2016 (Miyabuchi et al. 2018).During these periods, the crater lake had little to no lake water.It was restored in 2017-2018 and then depleted again in May 2019, which led to an ash eruption in 2019-2020 (Miyabuchi et al. 2021).Even after the cessation of the 2019-2020 eruption, the lake water did not return for more than a year.Subsequently, two phreatic eruptions occurred in October 2021 (Miyabuchi et al. 2022).In this study, we focused on the volcanic activity of the Nakadake crater from the end of the 2019-2020 eruption to the quiescent period after the 2021 eruptions.We divided the eruptive cycle into four phases, as described below (Fig. 2).

Phase 1: 2019-2020 eruption (July 2019 to early June 2020)
Ash eruptions began at the end of July 2019 at a main vent (the 191 vent) at the center of the crater floor (Fig. 1b).The ash eruptions continued until early June 2020 with occasional emissions of white plumes without ash (Miyabuchi et al. 2021).On June 8, 2020, an eruptive plume was emitted at a high temperature of > 300 ℃ (Fig. 3a).After a heavy rainfall of 200 mm/d on June 11, a white plume without ash was emitted at a low temperature of ~ 100 ℃ on June 16 (Fig. 3b), which corresponds to the end of the 2019-2020 eruption.
Phase 2: Inter-eruptive quiescent period (mid-June 2020 to June 17, 2021) In Phase 2, the crater had three sources of heat discharge: the 191-vent plume, fumaroles on the southern wall of the crater (SF), and a weak fumarolic field on the southern half of the crater floor (FF) (Fig. 1b).No lake water was present despite this being a quiescent period.The thermal activity in the crater remained quiescent during this phase, although the maximum temperature of the FF sometimes increased from approximately the boiling point to 200-300 ℃ in October 2020, January-February 2021, and early May 2021 (Japan Meteorological Agency 2021).The last increase was accompanied by an increase in the amplitude of the continuous tremors, and the Japan Meteorological Agency (JMA) alert level temporarily rose from one to two.However, our observations did not detect any major changes in the surface thermal activity.Fig. 2 A timeline of thermal activity inside the Nakadake crater.The forms of heat discharge (ash, plume, fumaroles) and their sources (191, BF, SF, NP, and CC) are plotted in each phase.The size of the plumes represents the relative magnitude of the heat discharge (not to scale).Annotations a-j indicate dates of shooting photographs shown in Fig. 3 in the NP on August 31 (JMA 2021).On September 1, the 191-vent plume became smaller than the NP plume.
On September 9, a UAS survey revealed that the water level in the NP rose to ~ 20 m below the crater floor.On September 22, hot water was present in the 191 vent as well as in the NP (JMA 2021).On October 4, the NP plume was still larger than the 191-vent plume (Fig. 3d), and the water level in the 191 vent rose further.From October 7, clear thermal manifestations in the crater began, although no our observation was conducted.On October 7, weak fountains occurred in the hot water pond in the 191 vent (JMA 2021).On October 12, there were escalated thermal manifestations over the entire crater; the water level of the 191 vent reached the elevation of the crater floor (20-m rise from September 9), the 191-vent plume visually increased in size, and the area of geothermal anomalies over the FF expanded (JMA 2021).From October 8, a strain meter at the Hondo Tunnel station began to observe radial compressional strain, which indicated inflation at a shallow depth beneath the crater.On October 13, the amplitude of the continuous tremors rapidly increased.
Phase 4: From eruptions to a quiescent state (October 14,

2021, to early March 2022)
A small phreatic eruption occurred at 4:43 JST (UTC + 9 h) on October 14, and another pit hole formed south of the 191 vent (hereafter referred to as the south pit hole: SP).Plumes of almost the same size were emitted from three sources: the 191 vent, the SP, and the NP (Fig. 3e).After October 15, no plume was observed in the SP (Fig. 3f ).The amplitude of continuous tremors began to increase again on October 18.Consequently, a more intense phreatic eruption occurred at 11:43 JST on October 20, with a plume height of ~ 3500 m.This eruption produced a low-temperature pyroclastic density current and ballistic ejection, with a maximum reaching distance of 900 m (Yokoo et al. 2022).This eruption caused the 191 vent and NP to merge into a single hole at the center of the Nakadake crater (CC).The crater floor deepened by ~ 40 m, and muddy water was present in the CC (Fig. 3g).
After the two eruptions, the CC continued to emit an enormous plume until early March 2022.The plume emissions from the SF intensified after the eruption (Fig. 3h) and occasionally merged with the CC plume.SO 2 gas flux also increased substantially compared with those in Phase 3 (300-800 t/d), reaching a maximum flux of > 5000 t/d on some days (JMA 2021).Because of the strong plume emissions (Fig. 3h and i), the inside of the crater could not be clearly observed on many days.In early March 2022, the strength of the plume emissions declined, and the SO 2 flux decreased to < 1000 t/d, which allowed the crater inside to be observed.The surface thermal activity returned to quiescence with the restoration of the crater lake (Fig. 3j), and the intensified plume emissions caused by the 2021 eruptions ended in early March 2022.

Thermal observation and analysis
The heat discharge sources in the Nakadake crater during Phases 1 to 4 were the plumes of the 191 vent, the NP, the SP, and the CC, and the fumarolic fields of the FF and SF (Fig. 1b).The HDRs of vent plumes ranges from 300 to 5000 MW, whereas the total of the FF and SF is up to ~ 20 MW (Nashimoto and Yokoo 2023).Thus, the heat discharge over the crater is dominated by plume-laden heat (> 93% of the total).To determine the magnitude of HDR in the crater, we focused on the 191 vent in Phases 1-2, the 191 vent and the NP in Phase 3, and CC in Phase 4.During Phases 1-3, we mainly observed the crater inside from the crater rim by using TIR cameras.During Phase 4, we mainly observed from outside the crater at a distance using visible cameras.

Thermal infrared (TIR) camera observations
We observed plume temperatures using TIR cameras from May 2020 to October 2021 (Phase 1-3).Three observation points were set up on the southwest (A), northwest (B), and west (W) sides of the crater rim (Fig. 1b).The observation was conducted once to twice per month, except in Phase 3. We used a handheld TIR camera of InfRec G120EX (Avionics Japan) to observe the plumes of the 191 vent and NP (Fig. 3a-f ) at 1/3 frame per second (fps) from May to December 2020 and at 10 fps from January 2021 to October 2021.We mainly made observations between sunset and dawn to avoid the influence of solar radiance.Because the recorded brightness temperature is affected by atmospheric absorption, object emissivity, and plume transparency (e.g., Harris 2013;Gaudin et al. 2016), we corrected for these factors to retrieve plausible plume temperatures.The atmospheric absorption was corrected by using an internal algorithm implemented in the G120EX.The emissivity of the plume with the water droplets was assumed as 0.97 (e.g., Gaudin et al. 2016).Plume transparency was assumed as zero because we focused on the opaque white part of the plume.
In addition to the handheld TIR cameras, we also performed aerial TIR observations on September 23, 2021 to focus on the plume and hot water in the NP.We used a TIR camera (FLIR Zenmuse XT2) mounted on a UAS (DJI Matrice 200).Because an XT2 image is in FLIR systems format, which allows for access to raw radiant energy records, an end-user can correct for the effects of atmospheric absorption and emissivity using an arbitrary method (Yokoo and Ishii 2021).The emissivity of the water surface and condensed plume was set to 0.97, just like for G120EX.The atmospheric absorption effect was calculated by using a LOWTRAN-based model (Minkina and Dudzik 2009).The observed temperatures of the water in the NP and of both plumes (the 191 vent and NP) were used to calculate the total HDR of in the crater.To calculate the HDR derived from the hot water in the NP, we considered realistic values of seepage flux and evaporative heat estimated from the water surface area and temperature (Additional file 1: Text S1; Table S1).
Because we used different cameras (G120EX and XT2) and correction methods for the TIR images, the measured temperatures differed even when the same target was captured.To estimate the temperature difference between the handheld and the aerial cameras, we took images of the crater lake as a substitute for the vent plume.As the emissivity of the target object (water) and the line-of-sight distance from the observation points to the crater lake surface (250-300 m) were almost the same as those of the plumes in the crater, we speculate that the correction obtained for the lake is approximately applicable to plume temperature.From several measurements in 2022, the XT2-based temperature was found to be 9 ± 2 ℃ lower than G120EX-based temperature (Additional file 1: Table S2).Both corrected temperatures were still lower than the true plume temperature because we assumed spatially homogeneous humidity.Thus, we corrected the XT2-based temperature by adding 9 ℃ first of all, and used this corrected plume temperature to calculate the plume HDR.

Estimation of HDR retrieved from TIR images
TIR cameras recorded the temperature of a plume, which is a mixture of volcanic gas emitted from a vent, and entrained ambient air.Given the temperatures of the plume and ambient air, the HDR and corresponding mass flux of volcanic gas included in the plume can be estimated by considering the energy and mass budgets of the plume, volcanic gas, and ambient air (Fukui 1995;Matsushima et al. 2003;Matsushima 2005;Witter et al. 2012;Gaudin et al. 2016).In this study, we applied the method proposed by Matsushima (2005) to estimate the HDR from TIR images.Volcanic gas is emitted from a vent, mixed with the entrained ambient air, and then cooled to a temperature recorded by the TIR cameras.Assuming volcanic gas consists of only H 2 O, the mass conservation law of water and air included in each thermal component (volcanic gas, ambient air, and plume) can be expressed as where m is the mass flux (kg/s) of each component and phase.The subscripts gas, air, and plu indicate volcanic gas, ambient air, and plume, respectively.The superscripts v, l, and a indicate water vapor, liquid droplets, and dry air, respectively.Under ideal gas approximation, m v plu , m v air , and m a plu can be expressed as (1) where e is the saturation vapor pressure, and P is the total pressure, which is assumed as the atmospheric pressure (0.89 × 10 5 Pa at the elevation of the crater floor).R is the gas constant ( R v =461.5 J/kgK, R a =287.04J/kgK), ω is the mixing ratio, T plu is the recorded plume temperature, RH is the relative humidity, v is the flow velocity, and r is the plume radius.The system of equations is closed by considering the conservation law of thermal energy: where h is the specific enthalpy of each component and phase.The left-hand terms represent the total enthalpy of the volcanic gas and ambient air before mixing, whereas the right-hand terms represent the enthalpy of the plume.In Eqs.(1-6), the remaining unknowns are h v gas and m v gas , and only h v gas m v gas (HDR) is uniquely determined.Because h v gas can be constrained by the maximum temperature of the vent exit observed by TIR cameras, m v gas (H 2 O flux) can be estimated separately.
The input parameters were obtained from TIR images and in-situ meteorological measurements.h v gas can be calculated by referring to the maximum temperature of the plume at a vent exit.Other specific enthalpies ( h v air , h a air , h v plu , h l plu , h a plu ) can be calculated from temperatures and specific heat of water vapor, liquid water, and dry air.The temperature and relative humidity in ambient air were obtained by in situ measurements at the crater rim.The relative humidity within the plumes was assumed 100%.We considered a profile perpendicular to the direction of the plume flow and defined the average temperature of this profile as T plu (Fig. 4a).The plume radius r was retrieved from the profile with temperature deviated from ambient temperature (Fig. 4b).To estimate the plume-rise velocity v , we employed an image velocimetry technique called optical flow.Several studies have demonstrated that optical flow accurately measures the velocity fields of SO 2 gas with high accuracy (e.g., Peters et al. 2015;Gliß et al. 2018) and is also applicable to TIR images (Lopez et al. 2014).We used the Farnebäck algorithm (Farnebäck 2003), a well-established optical flow algorithm, implemented using the OpenCV Python library.This algorithm has seven parameters for velocity analysis that need to be optimized for a target field to retrieve the velocity with sufficient accuracy.We used particle-image velocimetry standard images with known mean velocities of 2.5-22 pixels/frame (Okamoto et al. (4) RH air e v air (P−RH air e v air ) m a plu (6) to optimize parameter setting (Additional file 1: Table S3), which reproduced the mean velocity with > 98% accuracy as long as the mean velocity of the target field was < 10 pixels/frame.We applied this method to the TIR images acquired in 2021.We created inter-frame differences of consecutive TIR images and applied the Farnebäck method.We stacked the estimated velocity fields for time-averaging (Fig. 4c) and took the mean value of the profiles in the direction orthogonal to the flow direction as the mean velocity (Fig. 4d).For TIR images acquired in 2020, sampled at 3-s intervals, we estimated the velocity by manual tracking (Ilanko et al. 2019;Tamburello et al. 2019) because the displacements of plumes between frames were too large to apply the optical flow.

Plume-rise observation and analysis
We made both ground-based and aerial observations for the plume rise using visible video cameras.Ground-based observations were made from May 2020 to February 2022 at Kusasenri and Sakanashi stations (Fig. 1a) using cameras at 1-10 s intervals.An aerial observation was made on January 12, 2022, from 1 km southwest of the crater at 3-s intervals (Fig. 1a) using a camera (XT2) mounted on a UAS (M200).To retrieve the HDR of the plumes, we applied the Briggs model of buoyant plumerise (Briggs 1969).This model assumes that a bent-over plume rise is driven only by buoyancy under a constant crosswind and neutral atmospheric conditions.The HDR of a plume (J/s) is expressed as (Kagiyama 1981) where u is the wind velocity, C is a coefficient related to the plume flow shape, and α is a coefficient related to the meteorological condition and an entrainment coefficient.C is expressed as C = zx −2/3 where z is the altitude from a vent and x is the horizontal distance from the vent.C and u were directly estimated from plume images, and α was calculated as 2.4-2.7 × 10 4 based on in-situ meteorological parameters and a representative entrainment coefficient of 0.6 for a bent-over plume (Weil 1988).The HDR can be converted to H 2 O flux (kg/s) by dividing it by the specific enthalpy of water vapor at 96-200 ℃ (2667-2900 kJ/kg).The plotted estimation errors of the HDR (1σ) reflect the estimate uncertainties of wind direction, C , and u.
We estimated the optimal wind directions from meteorological data before estimations of u and C. We referred to anemometer records at Hondo station and JMA grid point values of the Meso-Scale Model at isobaric surfaces of 800-850 hPa, which corresponded to geopotential heights of 1500-2000 m.Because some wind data had (7) Q = αu 3 C 3 wind directions opposite to the actual direction predicted from plume images, we excluded these data from the direction determination and estimated the optimal wind directions from the remaining wind data.The winddirection correction yielded HDR values larger than those obtained without the correction (Fig. 5a), particularly on days when the wind blew into the depth direction of the image plane (Fig. 5b).The arial observation in January 2022 allowed us to observe the plume from a direction perpendicular to the actual wind direction (Additional file 1: Fig. S1); this greatly reduces the uncertainty of the wind direction and improved the estimation accuracy of the velocities and shape coefficients.Using a UAS yielded a small estimation error of the HDR on January 12, 2022; the ratio of the standard error to the mean was only ~ 30%, which was the smallest in Phase 4.
The shape coefficient C was estimated from consecutive images as follows.First, we extracted points on a plume streamline in a camera coordinate (Fig. 5c).Then, we projected the selected points onto planes along the wind directions in Cartesian coordinate (Fig. 5b) (Yanagisawa et al. 2022).Finally, we fitted a theoretical curve to the streamline points (Fig. 5c) and obtained the optimal value of C.
The wind velocity was estimated in a similar manner.First, we estimated the velocities within an xz-plane perpendicular to the camera's line-of-sight direction by using the zero-mean normalized cross-correlation of the brightness time series at two neighboring pixels on the plume streamline (Fig. 5c and d) (Crone et al. 2008).We then corrected the effect of the wind direction to project the line-of-sight-plane velocities onto the wind-direction planes.

Multi-GAS observation
To validate our estimates, we compared them with gas observation obtained by differential optical absorption spectroscopy (DOAS) and multi gas analyzer system (Multi-GAS) (Shinohara 2005).JMA conducted DOAS to measure SO 2 flux at least twice a month.The Multi-GAS instrument was installed in the western part of the crater rim (near station W) and continuously measured gas compositions, including H 2 O/SO 2 and CO 2 /SO 2 (Morita et al. 2022) during August 2019-August 2021 (corresponds to Phase 1 to 3).Because this Multi-GAS instrument stopped operating two months before the 2021 eruptions, we conducted UAS-mounted Multi-GAS observation on January 12, 2022, three months after the eruption (Phase 4), to measure gas composition of posteruptive plume. For

Heat discharge rate
Because we utilized two methods to estimate the plume HDR, we needed to confirm their consistency.We compared the HDRs with observation time intervals of < 12 days.This comparison showed no clear bias between them and that they were consistent with a deviation by a factor of 1/3 to 2 (Fig. 6).This consistency indicates that the HDR results can be plotted as a single time series to discuss the temporal changes.The estimation errors of the plume-rise analysis are consistently larger than that of the TIR-based modeling (Fig. 6) because the estimation error of the plume-rise parameters ( u and C ) propagates by the power of three (see Eq. 7).Thus, to discuss the temporal change in the plume HDR, we focused on the TIR-based results during Phase 1-3, with less scattered values and on the plume-rise-based results during Phase 4 because of the lack of TIR images observation.
Figure 7 shows the time series of the HDR and plume temperatures during 2020-2022 (see Table S4 for daily HDR).At the end of Phase 1 (May 13-June 8, 2020), a large HDR of 1-4 GW was observed (Fig. 7a), which may be plausible because of the vigorous emission of ash-laden plumes with temperature of > 300℃ (Figs.3a  and 7b).In Phase 2, the HDR decreased to 600-800 MW on June 16, 2020, which corresponds to a decrease in plume temperature to 100 ℃ (Fig. 7b).From June to December 2020, the HDR remained almost constant at 450-800 MW.From January 2021, it decreased and remained at 300-550 MW until mid-June 2021.In Phase 3, the HDR of the 191-vent plume was slightly lower (~ 250 MW) on June 23, 2021 (Phase 3), five days after the formation of the NP.Because the NP had already emitted a fumarolic plume of almost the same size as the 191-vent plume on this day (Fig. 3c), the total HDR of the crater could be approximately twice the HDR of only the 191-vent plume (~ 500 MW).
On September 23, 2021, the last observation before the 2021 eruption, the HDR was 260-350 MW at the NP and 70-160 MW at the 191 vent; its total was 330-510 MW (Fig. 7a and b).Considering the uncertainty for the corrected temperature of the two TIR cameras (7-11 ℃), the total HDR had a range of 400-600 MW (Additional file 1: Fig. S2).These values are consistent with those before the depression (300-550 MW) and just after the depression (~ 500 MW).The conservation of the total HDR before and after the depression implies that the 191 vent and the NP were connected to a common conduit at depth.The maximum temperature of the fumarolic plumes was ~ 90 ℃, which had remained constant since January 2021 (Fig. 7b).On this date, the NP hosted both hot water and a fumarolic plume, which suggests that volcanic gas supplied from the pit floor maintained the hot water at a temperature of ~ 90 ℃ near boiling point of 96 ℃.The heat input was estimated to be 30-40 MW, which was < 10% of the total HDR in the NP (Additional file 1: Text S1).Thus, there was no evident change in the HDR up to three weeks before the October 14 eruption.As a long-term trend during Phases 2-3, the HDRs were consistently higher than those during past quiescent periods with a crater lake (150-260 MW) (gray-shaded zone in Fig. 7a).
On October 14 (5 h after the first eruption), the HDR was ~ 800 MW from the plume-rise analysis (Fig. 7c).Although plumes were emitted from three sources, the 191 vent, SP, and NP (Fig. 3e), we could not identify the source of the plume in the analyzed images.It is possible that we missed other plumes that were not captured in the analyzed images, which would mean that the total HDR was underestimated.On October 15, while plume sources were both the 191 vent and NP, we could only estimate the HDR of the 191-vent plume with the TIRbased method (540-650 MW).The total HDR may be up to twice this value at 1100-1300 MW (Fig. 7c) because the two plumes were almost the same size.Thus, the HDR increased by 2-4 times from September 23 (Fig. 7c).Hereafter, only plume-rise analysis was available for estimating the plume HDR.On October 19, the HDR further increased to 2.6-5.3GW, the largest value since the end of the 2020 eruption, which coincided with an increase in SO 2 flux (Fig. 7c).Around 16:00 JST on October 20, 5 h after the second eruption, a significant HDR was also estimated as 1.7-4.4GW (Fig. 7c).It decreased to 1.1-2.4GW on October 28 and remained almost constant at around 1 GW for the rest of Phase 4, except for November 28, 2021 and February 24, 2022.
On both days, the HDR values temporarily increased in a spike-like manner (Fig. 7a).On the UAS-based observation on January 12, 2022, the HDR (0.8-1.6 GW) is the best constrained during Phase 4.

Mass flux
Figure 8a shows the temporal changes in the H 2 O flux of volcanic gas included in the emitting plumes.In Phase 1, the H 2 O flux was estimated as 60,000-110,000 t/d.In Phase 2, it decreased from 15,000-25,000 t/d in the first half (June to December 2020) to 10,000-17,000 t/d in the second half (January to June 2021).This is consistently higher than that emitting during the past quiescent periods of 2000-2003and 2006-2009 (3500-11,000 t/d;Saito et al. 2008;Terada et al. 2012).In Phases 2-4, the temporal characteristics of H 2 O flux was almost the same as those of the HDR, which also decreased from Phase 1 to 3 and increased from Phase 3 to 4 (Fig. 7).This is because the maximum plume temperature was almost constant, at around 100 ℃, resulting in a temporally constant gas enthalpy of ~ 2700 kJ/kg.The time series of H 2 O and SO 2 fluxes showed roughly the same temporal variation, except in Phase 1 (Fig. 8a), which implies that H 2 O/SO 2 remained almost constant throughout Phases 2 to 4. This correlation was also found before the phreatomagmatic eruption in 2016 (Morita 2019).The H 2 O/SO 2 may be systematically higher in Phase 1 than in Phases 2-4, which indicates that the gas composition during the 2020 ash eruption was unrealistically poor in SO 2 compared to the composition during the inter-eruptive period (Phases 2-3).

Validation of the H 2 O flux through comparison with gas observations
We validated the analyzed results of our thermal observations by comparison with gas observations (DOAS and Multi-GAS) in terms of H 2 O/SO 2 and H 2 O flux.In Phase 1, the H 2 O/SO 2 was 10-90 for Multi-GAS and 100-400 for the thermal and DOAS observations (Fig. 8b).The former values were measured during ash eruptions and were consistent with those observed at the 2015 ash eruptions (30-60) (Shinohara et al. 2018).In contrast, the latter values were much larger than the Multi-GAS measurements and also even 2-4 times larger than those for inter-eruptive periods (Phases 2 and 3).This indicates that the H 2 O/SO 2 calculated by the thermal and DOAS observations in Phase 1 was counterintuitively poor in SO 2 despite Phase 1 being an eruptive period.This higher H 2 O/SO 2 can be attributed to either overestimating the H 2 O flux or underestimating the SO 2 flux.
The possibility of the overestimation of H 2 O flux can be rejected as follows.In Phase 1, the plume temperature was higher than 300 ℃, which resulted in the larger TIRbased HDR in Phase 1 than the HDR in Phases 2 and 3 at a lower temperature of ~ 100 ℃ (Fig. 7b).Furthermore, the temperature recorded with TIR cameras may have been underestimated because the emissivity of the ashy plumes was not be fully corrected.Thus, these did not result in an overestimation of the H 2 O flux obtained from the TIR images.On the other hand, the plume-rise analysis should have underestimated the H 2 O flux because we neglected the thermal contribution of ash included in the plume.However, this contribution was only ~ 2 MW when assuming a specific heat of ash of 1 kJ/kg K, a temperature difference of 400 K, and ash mass flux of 500 t/d (Miyabuchi personal communication 2020).This HDR is negligible compared to the HDR magnitude of ~ 1 GW.Furthermore, the HDR and H 2 O flux estimated by both methods were consistent in Phase 1 (Figs.6-8).Thus, the H 2 O flux based on our thermal observations in Phase 1 is unlikely to be overestimated.
Instead, the underestimation of SO 2 flux is likely because the presence of ash in eruptive plumes can decrease the amount of ultraviolet rays absorbed by SO 2 gas (Andres and Schmid 2001;Kazahaya et al. 2013).In addition, SO 2 flux frequently increased by a factor of two to three in Phase 1 (Fig. 8a).However, such rapid increases could not be reproduced by a linear time interpolation of the SO 2 flux (see the section 'Multi-GAS observation'), and the interpolated SO 2 flux on the dates of thermal observations was likely to be frequently underestimated by half to one-third of the actual flux.Therefore, the unrealistically high H 2 O/SO 2 in Phase 1 according to thermal and DOAS observations was probably due to underestimate of the SO 2 flux.
In Phases 2 and 3, the Multi-GAS measurements showed H 2 O/SO 2 of 50-150, systematically higher than those in Phase 1 (Fig. 8b).This was probably due to a decrease in the contribution of magmatic gases and a shift to a gas composition influenced by a hydrothermal system.The H 2 O/SO 2 calculated from the thermal and DOAS observations varied within 30-300 with an average of ~ 100, which was consistent with Multi-GAS measurements.
In Phase 4, H 2 O/SO 2 measured by the UAS-mounted Multi-GAS was 58, whereas the H 2 O/SO 2 calculated from thermal and DOAS observations varied within 20-200, with an average of ~ 70 (Fig. 8b).These results suggest that H 2 O/SO 2 decreased from Phases 2-3 to Phase 4, which corresponds to an increase in the contribution of magmatic gas.This is also supported by the large SO 2 flux in Phase 4 (2000-6000 t/d), which was comparable to that of the 2019-2020 ash eruption.Although comparison between the Multi-GAS and thermal monitoring results in Phase 4 is possible for only one epoch (January 12, 2022), its validity is plausible because all observations (visible camera, DOAS, and UAS-based Multi-GAS) were successfully conducted on this date, which allowed for comparison without any time interpolation.On this date, the H 2 O fluxes and H 2 O/SO 2 estimated by thermal and gas observations showed good consistency; H 2 O flux was 28,000-55,000 t/d for UASbased plume-rise analysis and 24,000-40,000 t/d for DOAS and UAS-based Multi-GAS (Fig. 8c); meanwhile H 2 O/SO 2 was 42-104 (average 73) for the plume-rise analysis and DOAS and was 58 for Multi-GAS (Fig. 8b).While H 2 O/SO 2 was almost constant in Phase 4, it temporarily increased on November 28, 2021, which corresponded to a pulse-like increase in H 2 O flux (Fig. 8).This increase may be due to a discrepancy between the timings of the thermal and DOAS observations.Because the plume emissions in Phase 4 were unstable and increased intermittently, interpolated SO 2 flux on dates of the thermal observation could not fully reproduce true temporal change, as mentioned in Phase 1.
The continuous Multi-GAS station might have sampled both the 191 vent and SF plumes.However, the good consistency between the gas and thermal observations that mainly focused on the 191-vent plume implies that the primary gas source was the 191 vent.Furthermore, UAS-based measurements of each gas source on October 20, 2020, showed that CO 2 /SO 2 was 2.0 near the 191 vent and 0.8 at the SF (Tsunogai et al. 2022).The CO 2 /SO 2 of the 191-vent gas (2.0) was close to the Multi-GAS values of 2.1-4.0 for August-December 2020 (Fig. 8d), which also supports the possibility that the Multi-GAS reflects the 191-vent gas composition.

Thermal activity and its association with the two eruptions
Below, we interpret the temporal evolution of heat discharge in the crater and its association with the two phreatic eruptions in 2021 by comparison with other geophysical and geochemical observations in each phase.

Phase 1-2: Decline of heat discharge at the end of the 2020 eruption
During the transition from Phase 1 to 2 (June 2020), the HDR dropped sharply from 2.5 to 0.7 GW (Fig. 7a), which coincided with SO 2 flux dropping from 2000 to 500 t/d (Fig. 8a).In addition, the total geomagnetic force changes measured around the crater reversed from demagnetization to magnetization around April-May 2020, indicating cooling at shallow depth (Kyoto University 2021).This coincides with the timing of the decrease in the surface heat discharge, which suggests a decrease in the heat supply from a greater depth, probably a magmatic region, into a shallow part beneath the crater.

Phase 2: Decreasing trend of HDR
From the first half (May 2020 to December 2020) to the second half (January 2021 to June 2021) of Phase 2, the average HDR decreased from 600 to 400 MW, which is coincident with a decrease trend of SO 2 flux (Fig. 9a).In early May 2021, the SO 2 flux and continuous tremor amplitudes increased intermittently (blue-shaded regions in Fig. 9a), which corresponded to an unrest due to enhanced degassing.In contrast, no such increase was observed for HDR, which can be likely attributed to the infrequency of the thermal observation; the enhanced SO 2 gas emissions on May 11 and 17 were observed between the timings of thermal observations on May 6 and 19.Therefore, the thermal observation could not capture such short-period unrest.H 2 O/SO 2 showed a sharp decrease at the unrest (Fig. 8b), but this was also due to a similar reason; the interpolated SO 2 flux on the day of our observation was high at 1000 t/d.Because no DOAS was conducted on this day, the validity of the interpolated SO 2 flux around the intermittent SO 2 flux increase is ambiguous, and we cannot confirm that such a low H 2 O/SO 2 is reliable.

Phase 2: Origin of water in the plumes
We here assess how much magmatic gases contribute to the observed plume through mass and energy budget calculation.In Phases 2 and 3, the plume temperatures of the 191 vent and the NP were low at < 120 ℃, close to the boiling point (Fig. 7b).On the other hand, the apparent equilibrium temperature of gas sampled near the crater floor was ~ 630 ℃ in October 2020, which was interpreted as the temperature before cooling by near-surface cold water (Tsunogai et al. 2022).These observations imply that groundwater or hydrothermal water effectively cooled the higher temperature magmatic gases.This situation is expressed by mass ( m ) and energy ( hm ) conservation per unit time as where the subscripts w and g correspond to liquid water and magmatic gases, respectively.From Eqs. ( 8) and ( 9), the mass ratio of magmatic gases to the plume m g /m plu is expressed as where h plu was 2677-2900 kJ/kg in Phase 2, h g is assumed to be 3700 kJ/kg, corresponding to water steam at 630 ℃ (Tsunogai et al. 2022), h w was assumed as 420 kJ/ kg, corresponding to liquid water at 100 ℃, and m plu was estimated as 170-290 kg/s in the first half of Phase 2 (Fig. 8a).
We focus on energy and mass budget in the first half of Phase 2 (June to December 2020) because there is a reasonable estimation of apparent equilibrium temperature of magmatic gas in October 2020.By using Eq. ( 10), we obtained m g /m plu of 0.69-0.71,which indicates a large portion of plumes was derived from magmatic gases.The estimated magmatic gas flux (120-205 kg/s) was more than twice that estimated in 2006-2009 (40-65 kg/s) (Terada et al. 2012).This shows that magma degassing was more active in Phase 2 than in the past quiescent period, even though the volcanic activity in Phase 2 was seemingly quiescent, such as low seismicity and no ground inflation (Fig. 9a).In addition, (8) the CO 2 flux was 700-4000 t/d in Phase 2 (Fig. 8e), larger than that of the past quiescent period (502-692 t/d; Saito et al. 2007), which suggests that degassing from a deep magmatic region (> 10 km) was enhanced compared to degassing during the past quiescent periods.
This degassing requires a substantial amount of magma that can be preliminarily calculated from the magmatic H 2 O flux and H 2 O content estimated by melt inclusion analysis.We used the magmatic H 2 O flux of 120-205 kg/s and assumed the H 2 O contents of 0.5-2.26wt% estimated from eruptive material of the 2014-2015 eruptions (Saito et al. 2018;Kawaguchi et al. 2021).The magma mass required for the observed H 2 O gas discharge was estimated as 5200-41,000 kg/s, which is comparable to the estimated mass of 13,000-54,000 kg/s from melt inclusion analysis of sulfur content in the erupted scoria from the 2014 eruption (Saito et al. 2018).This consistency supports, in terms of thermal observation, the existence of density-driven magma convection within the magma plumbing system of Aso volcano, as suggested by previous studies (e.g., Saito et al. 2018).The estimated magma flux corresponds to a volume flux of 2-16 m 3 /s with a magma density of 2500 kg/m 3 , which is also consistent with that of other persistently degassing volcanoes, such as Satsuma-Iwojima at 7.5-13 m 3 /s (Kazahaya et al. 2002) and Ambrym at 25 m 3 /s (Allard et al. 2014).Although our estimates are based on the H 2 O flux and H 2 O content, these estimates are consistent with the estimated magma fluxes based on the SO 2 gas flux and the sulfur content at other volcanoes.As for SO 2 and CO 2 , the same calculation can be conducted using the gas fluxes observed in Phase 2 and the volatile content of the pre-degassing magma.Based on the SO 2 flux during Phase 2 (500-1000 t/d) and the sulfur content estimated in the 2014-2015 eruption (0.02-0.07 wt%; Saito et al. 2018), the pre-degassing magma flux can be calculated as 4100-29,000 kg/s.Based on the CO 2 flux during Phase 2 (700-4000 t/d) and the CO 2 content estimated in the 2014-2015 eruption (0.04-0.5 wt%; Saito et al. 2018), the pre-degassing magma flux can be calculated as 1600-120,000 kg/s.These show that the estimated magma fluxes are roughly in agreement with each other.Therefore, at least in Phase 2, the enhanced magma degassing continued to be as active as in the past eruption, even after the cessation of the 2019-2020 eruption.

Phase 2-3: Cause of the crater lake depletion
The primary form of heat discharge in the Nakadake crater differed depending on whether a crater lake was present; plume emission when the crater lake was absent (2020-2022) and evaporation from the lake surface when the crater lake was present (typical quiescent period).During the quiescent periods, the product of the specific enthalpy of water vapor and the volcanic gas flux emitting from subaqueous vents ( h v gas m v gas ) was estimated (Saito et al. 2008;Terada et al. 2012); this is the same as the HDR defined in our study, which allowed for a direct comparison of them.
HDRs during Phases 2 and 3 were constantly larger than during 2000-2009 (gray-shaded zone in Fig. 7a).Mass and energy budget modeling of the lake have shown that the presence of lake water is controlled by a balance between the inflow of volcanic gas and meteoric water and the outflow of evaporation and seepage (Terada and Hashimoto 2017).Based on the crater topography and hydrological settings, an addition of heat to the background thermal activity can accelerate water loss by enhancing evaporation from the lake surface, which depletes the lake water (Terada and Hashimoto 2017).Thus, the depletion of the lake water in 2020-2021 may be attributed to unstable conditions; a shallow part beneath the crater was overheated by heat addition (300-800 MW) anomalously larger than that for the quiescent periods (150-260 MW).The large heat supply for 2020-2021 was likely due to enhanced magma degassing of a shallow magma reservoir.This is suggested by contraction of the Choyo-Hondo GNSS baseline from early 2020 to August 2021 (Fig. 9), which corresponds to the deflation of a magma reservoir at a depth of ~ 5 km (Sudo et al. 2006).The reservoir deflation could be caused by a discharge of excess volatiles during magma convection within the conduit, which would have contributed to overheating the crater and hindered the restoration of the lake water.

Phase 3: Run-up process of the 2021 eruption
In Phase 3, the results of the thermal and gas observations contrasted with other geophysical signals, including seismicity, geomagnetic changes, and ground deformation.From June to the end of September 2021, no clear increase in the HDR was recognized, which is consistent with the lack of an apparent increase in the SO 2 flux (Fig. 9a).On the other hand, contraction of two GNSS baselines (C-H and H-S baselines) that had continued from the end of Phase 1 stopped, and they began to extend around August-September 2021 (Fig. 9b).These baseline extensions corresponded to the inflation of the magma reservoir at ~ 5-km depth (the C-H baseline) and a shallow hydrothermal system (the H-S baseline).Although the onset time for the extension of the H-S baseline was somewhat unclear because of the low signal-to-noise ratio, the baseline reversed from contraction to extension around September 2021.Furthermore, demagnetization began around August (Kyoto University 2021), which indicated heat accumulation at a shallow depth beneath the crater, and the continuous tremor amplitude began to increase in September to October (Fig. 9a).These geophysical signals consistently suggest an enhanced heat supply from a magmatic region deeper than the magma chamber at 5-km depth to the shallow hydrothermal system, in contrast to the constant thermal activity and SO 2 gas emissions.
Such behavior can be interpreted as partial sealing of a volcanic conduit (e.g., de Moor et al. 2016;Mick et al. 2021), which would lead to pressurization and heating beneath a seal, while suppressing thermal manifestation at the surface (e.g., Christenson et al. 2010;Tanaka et al. 2018).Given the subsurface structure beneath the Nakadake crater, this scenario likely explains the run-up process of the 2021 eruption.Based on the resistivity structure, a low permeability zone filled with sulfides and sulfur minerals was present from the surface to a depth of ~ 400 m, with a fracture-rich conduit zone (Kanda et al. 2019).In this situation, a temperaturedependent increase in viscosity of sulfur and/or the precipitation of sulfide and sulfur minerals is essential for seal formation (Hurst et al. 1991;de Moor et al. 2019).This idea is supported by the abundance of sulfur and hydrothermally altered minerals in the eruptive material of the 2021 eruption (Miyabuchi et al. 2022).In addition, sources of demagnetization and pressurization in past eruptions were emplaced ~ 200 m below the crater (Tanaka 1993;Ishii et al. 2023), and these sources are likely to be responsible for signals observed before the 2021 eruption.

Phase 4: Cause of the 2021 eruption
In contrast to the steady HDR until the end of September 2021, the thermal activity over the crater floor increased for October 7-12.This was accompanied by a strain change at Hondo station, indicating the inflation at a shallow depth.This thermal enhancement with shallow pressurization possibly suggests leakage of gas that had accumulated beneath the seal due to a further increase in heat supply from a greater depth, as reported for the Augustine 2006 eruption (Zhan et al. 2022) and the Poás 2017 eruption (de Moor et al. 2019).Consequently, on October 14, the onset of the first eruption in 2021 may have been caused by a rupture of the hydrothermal system, which was triggered by a further gas supply into to the already pressurized hydrothermal fluid beneath the seal.This sequence corresponds to the phreato-vulcanian eruption (Stix and de Moor 2018).
While the first eruption could have been caused by conduit sealing and further heat supply, the second eruption may have been caused by a substantial injection of magmatic fluid into the hydrothermal system, based on continuous pressurization of the magma chamber without relaxation even after the first eruption.As evidence, there was no clear co-eruptive deflation during the first eruption, and both the deep and shallow parts continued to be pressurized (Fig. 9a).In addition, the HDR and SO 2 flux increased from October 15 to 19 and reached their highest values since the end of the 2020 eruption (Fig. 7c).These imply that a massive gas continued to be supplied from the magma into a shallow part of the hydrothermal system even after the end of the first eruption.This further pressurized the stillpressurized hydrothermal system, leading to the larger eruption on October 20.
This sequence is similar to that of the phreatomagmatic eruption in October 2016, although the eruption intensities differed.Three months before the 2016 eruption, the magma chamber began to inflate as suggested by the extension of the C-H baseline.Simultaneously, SO 2 flux increased responding to the volatile-rich magma supply (Morita 2019).It rapidly decreased and remained low at 500-1500 t/d, possibly due to a conduit sealing, but it reached 15,000 t/d on October 7, which was the largest ever observed at Aso volcano (Morita 2019).Finally, 4 h after a minor eruption, an enormous phreatomagmatic eruption occurred on October 8, with a plume height of ~ 11 km (Ishii 2018).Thus, the sequences of the 2016 and 2021 eruptions are similar in terms of degassing activity, ground deformation, and a preceding minor eruption.The 2021 eruptions may be a smaller version of the 2016 eruption.The 2021 eruptions were previously regarded as one of the eruptions that occurred under groundwaterrich conditions while the magmatic activity declined (Miyabuchi et al. 2022).However, regarding pre-eruptive thermal activity and ground deformation, the 2021 eruptions were not a by-product eruption that usually occurs during a decline in magmatic activity at the Nakadake crater.Instead, the 2021 eruptions were driven by a new magma supply event independent of the magma supply that drove the 2019-2020 eruption.

Phase 4: Post-eruptive relaxation toward quiescence
In Phase 4, a significant HDR of 0.4-4 GW continued four months after the second eruption (Fig. 9a).The SO 2 flux also remained at high flux of 1000-6000 t/d, which was comparable to the flux during the 2019-2020 ash eruption (Fig. 9a).In addition, CO 2 /SO 2 decreased to 1.7, which is nearly the same as the values during the 2019-2020 ash eruption (Fig. 8d) and implies enhanced contribution of shallow magma degassing (e.g., Battaglia et al. 2019;de Moor et al. 2019).These results suggest that the contribution of magmatic gas was higher than in Phase 2-3.This enhanced release of magmatic gas can also be supported by volatile-rich magma supply to the shallow magma reservoir and possibly by subsequent accelerated magma convection within conduit.The magma supply continued until December 2021 as indicated by the inflation of the magma chamber (Fig. 9).Post-eruptive thermal activity tends to fluctuate unstably, and temporary increases were observed on November 28, 2021, and February 24, 2022 (Fig. 9).Because these dates correspond to intermittent increases in the SO 2 flux and tremor amplitude, respectively, these changes may be plausible signals of temporary enhancement of heat discharge.
After February 24, 2022, the plume-rise analysis was no longer applicable because the assumptions of the Briggs model broke down.This model assumes that plume trajectories are governed only by buoyancy and homogeneous crosswinds.However, the plume trajectory after the end of February 2022 was controlled by the complex wind distribution in the crater, and the plume often dragged on the surface because the plume was thin with weak buoyancy.Therefore, the plume-laden HDR after March 2022 was unknown.We speculate that the plume HDR was < 300 MW because the lower limit of our analyzed period was ~ 300 MW.This speculation is consistent with the HDR of ~ 100 MW of volcanic gas injected in the lake bottom in March 2022, estimated from thermal modeling of the crater lake (Nashimoto and Yokoo 2023).
In Phase 4, we could not distinguish whether the plume origin was the CC and/or SF.We inferred that the plume taken in Phase 4 was basically a larger plume coming from the CC.However, because the plume emitted from the SF was also enhanced after the 2021 eruptions (Fig. 3h), the actual HDR may be the sum of the CC and the SF plumes.Therefore, if we miss the SF plume independently emitted on some days, the estimated HDR in Phase 4 may be underestimated and correspond to the lower limit of the total values.However, the total HDR should be less than twice the present estimation because the SF plume was constantly smaller than the CC plume.

Comparison with seismicity of long-period tremors
Previous studies have often discussed the degassing activity at the Nakadake crater in terms of its relationship with very long period seismicity (VLP), of which the local name is long period tremors (Sassa 1935;Yamamoto et al. 1999).VLP events are regularly observed regardless of eruption occurrence or cessation in Aso volcano (e.g., Kaneshima et al. 1996;Kawakatsu et al. 2000).Typical VLPs are mainly excited by deflation of a tensile crack source at a shallow depth of 0.3-2.8km beneath the crater (Yamamoto et al. 1999), and may result from fluid discharge from the crack (e.g., Kawakatsu et al. 2000;Ishii et al. 2019).The time series of the daily sum of the VLP's root-mean-squared (RMS) amplitude are plotted (Fig. 9a), which has a period of ~ 15 s.Because the amplitude of VLP is proportional to its moment (e.g., Maeda et al. 2019), the time series of the RMS can be regarded as a time function of the daily moment release of the VLP source, although the absolute value itself has no physical meaning here.The daily RMS decreased from Phases 1 to 3, then increased just before Phase 4, and finally decreased at the end of Phase 4 (Fig. 9a).These temporal characteristics are roughly consistent with those of HDR and SO 2 flux, which supports the interpretation of the previous study that VLPs are related to fluid discharge from the crack source (Kawakatsu et al. 2000).However, a closer look shows that these time series are not well correlated over the entire study period.In particular, the RMS drastically decreased from Phase 2 to 3, but no such decrease was observed in the corresponding HDR and SO 2 flux (Fig. 9a).In addition, the RMS in Phase 4 returned to the same level as in Phase 2, but the HDR and SO 2 flux were higher in Phase 4 than in Phase 2. These facts perhaps imply some correlation between the VLP amplitude and HDR and gas emissions throughout the study period, but the correlation coefficient changed with time (Fig. 10).In Phases 2 and 3, the RMS of the VLP events had a clear correlation with gas fluxes (H 2 O and SO 2 ).From Phase 2 to 3, the RMS decreased while the gas fluxes were nearly constant, which suggests some link between the crater depression and subsurface VLP excitation process, although the detailed mechanism is unclear here.In Phase 4, the correlation was still present, but it was weaker than in Phase 2-3.This is likely due to the transition of conduit system from relatively closed to open state before and after the eruption.The correlation change can also be attributed to other factors, such as changes in the conduit wall rock strength (Niu and Song 2020) and fluid flow path, which could switch degassing with or without VLP (Ichimura et al. 2018).The complex physical relationship between the shallow conduit system represented by the VLPs and the surface fluid discharge process, as captured by the thermal and gas observations, requires further analysis.

Conclusions
We evaluated the heat discharge rate of volcanic plumes associated with the 2021 eruptions of Aso volcano.Despite using two simple methods, our estimates were consistent with an independent dataset of volcanic gas monitoring.Based on our results, we identified transient thermal activity in each phase: Phase 1: A massive heat discharge with several GW was estimated at the end of the 2020 ash eruption.Phase 2: During the quiescent period between the 2020 and 2021 eruptions, the heat discharge rate decreased from in Phase 1 but remained at 300-800 MW, which was higher than that of the typical quiescent periods.This anomalously large heat discharge continued to overheat the crater floor for more than a year, probably hindering the restoration of Yudamari crater lake, which is typically present during quiescent periods.Phase 3: In the run-up stage, a comparison with steady heat discharge (~ 500 MW), ground deformation and demagnetization suggested that hydrothermal sealing caused subsurface pressurization and heating 1-2 months before the eruptions, and resulted in the first eruption.Phase 4: Further heat supply into the still-pressurized hydrothermal system after the first eruption, eventually resulted in the second eruption.During the post-eruptive period, a massive heat discharge of more than 1 GW continued for 4 months, which was supported by a continuous magma supply from the deeper magmatic region.At the end of Phase 4, the thermal activity in the crater returned to quiescence with the restoration of the lake water.Thus, we revealed the temporal evolution of thermal activity throughout a cycle of eruptive activity at Aso volcano.
We found three important aspects of thermal activity at Aso volcano; steady heat discharge during preeruptive quiescence, a precursory increase in heat discharge observed on the day before the second eruption, and similarity between the 2021 eruption and the 2016 eruption in terms of thermal, gas, and geodetic monitoring.However, we have not been able to quantify a short-term (~ 1 weak) precursory change in heat discharge prior to the first eruption.Capturing such short-term thermal manifestations with high temporal sampling is crucial for disaster mitigation of future eruptions with a short-time thermal precursor, as demonstrated by volcanic gas monitoring of the Poás 2017 eruption.

Fig. 1 a
Fig. 1 a Map of observation stations in Aso caldera.The circles indicate the locations of visible cameras (Sk: Sakanashi, K: Kusasenri, U: UAS hovering position on January 12, 2022).The cross indicates Sunasenri seismic station.The diamonds indicate GNSS stations (C: Choyo GEONET, H: Hondo, Sn: Sunasenri).The dashed lines (C-H and H-S) correspond to GNSS baselines shown in Fig. 9.The inset square corresponds to the region shown in (b).b Map of the Nakadake crater.The squares (A, B, and W) indicate the observation sites of the handheld TIR camera.The black circles indicate the hovering position of a UAS on September 23, 2021.The red circles indicate the vent formed in 2019 (the 191 vent), the north pit hole (NP), the south pit hole (SP), and the south wall fumaroles (SF).The red dashed line indicates the crater floor fumaroles (FF).The DEM was generated by merging the UAS survey result on September 9, 2021 with a GSI 10-m mesh DEM.The contours indicate elevation intervals of 3 m

Phase 3 :
From depression to the run-up of the 2021 eruptions (June 18, 2021, to October 13, 2021) On June 18, 2021, the crater floor was depressed and a pit hole was formed to the north of the 191 vent (hereafter, the North pit hole; NP).No fumarole or water was present in the NP on June 20, but a fumarolic plume was emitted on June 21.The plume grew to the same size as the 191-vent plume on June 23 (Fig.3c), and temporarily disappeared on July 27.The NP floor was visible on this day, and the depth was ~ 50 m, according to an unmanned aerial system (UAS) observation.Hot water was present

Fig. 3
Fig. 3 Photos of the thermal activity inside the Nakadake crater during Phases 1-4.The parentheses indicate photo shooting sites.a High-temperature (> 300℃) ashy plume emitted from the 191 vent on June 8, 2020.b Low-temperature (90℃) white plume emitted on June 16, 2020.c Inside the crater on June 23, 2021, which was 5 days after the depression.d Aerial view of the 191 vent and the NP on October 4, 2021, which was 10 days before the first eruption in 2021.e Inside the crater 5 h after the October 14, 2021 eruption.f Inside the crater on October 15, 2021.No plume was emitted from the SP.g View of the crater after the eruptions.A larger plume was emitted from the CC (December 10, 2021).h Plumes separately emitted from the CC and SF (December 10, 2021).i The UAS-based observation for plume-rise analysis (January 12, 2022).j Restoration of the crater lake representing a typical quiescent period (March 4, 2022)

Fig. 4
Fig. 4 Retrieval of plume parameters from TIR images taken at station W at 18:30 JST on March 23, 2021.a Stack of TIR images (30-s average) of the 191-vent plume.b Plume diameter and mean plume temperature T plu are obtained along a temperature profile.c Stack of velocity fields (30-s average) of the 191-vent plume.d Mean and standard deviation of the velocity are retrieved along ten profiles with the same length and distributed perpendicular to the flow direction (the black line in c) comparison with the H 2 O flux estimated from the thermal observations, we calculated the H 2 O flux from the gas observations by multiplying the H 2 O/SO 2 measured by Multi-GAS and SO 2 flux measured by DOAS.Because the timings of the DOAS and Multi-GAS observation do not always coincide, we calculated the SO 2 flux on the date of the Multi-GAS observations by linear interpolation of the original time series of SO 2 flux.Similarly, we calculated H 2 O/SO 2 from the H 2 O flux from thermal observations and the SO 2 flux and compared it with the H 2 O/SO 2 measured by Multi-GAS.The SO 2 flux on the date of thermal observations was interpolated in the same manner.The error values in the time-interpolated SO 2 flux were assumed as the temporal average of the ratios of the error values to the mean fluxes from 2019 to 2022, which was ~ 28%.Similarly, we calculated the CO 2 flux from CO 2 /SO 2 and the interpolated SO 2 flux.

Fig. 5
Fig. 5 Example of plume-rise analysis on February 18, 2022.a Distribution of the estimated HDR depending on wind directions (WD), measured in terms of angle clockwise from north.The optimal HDR values vary depending on whether the WD correction was included.b Positional relationships of the camera, vent, camera line-of-sight (LOS) direction, camera image plane, and WD.c 500-s stack of consecutive images of plumes.The black circles are pixels used for C curve fitting.The red line is the estimated C curve.p1 and p2 correspond to points shown in time series of pixel brightness in (d).d Velocity estimation using the zero-mean normalized cross correlation (ZNCC) of the time series of pixel brightness of neighboring two points (p1 and p2) along the time-averaged plume streamline as shown in (c)

Fig. 6
Fig. 6 Comparison between HDR values estimated from plume-rise and TIR-based model.Time dates for plotted TIR-based HDR were selected to be nearest to the dates of the plume-rise observation (time differences were less than 12 days)

Fig. 7 Fig. 8
Fig. 7 Results of thermal data analysis.a plume HDR based on TIR images (black circles) and visible images (white triangles).The gray-shaded zones indicate HDR range during the past quiescent period (Saito et al. 2008; Terada et al. 2012).The parentheses indicate the corresponding photo index (a-j) in Fig. 3.The blue arrow indicates the zoom-up period in (c).The smaller transparent circles and larger solid circles indicate individual estimates by the TIR-based method and their weighted average on each day. 1 The plume included a small amount of ash of < 500 t/d (Miyabuchi pers.com.) until June 8, 2020. 2 A strong wind blew parallel to the camera's LOS direction.This configuration could greatly underestimate a plume radius, which would result in an underestimation of HDR. 3 The plotted HDR value is twice the HDR of the 191-vent plume, which reflects the contribution of the NP plume. 4Plotted HDRs are the sum of the 191-vent and the NP plume. 5The plotted HDR is twice the HDR of only 191-vent plume, reflecting the contribution of the NP plume.b Time series of the maximum temperature of the 191 vent.The symbols in black, red, and blue indicate observation sites of A, B, and W, respectively.c HDR time series for September-October 2021.The red lines indicate the timings of the eruptions.Blue symbols indicate SO 2 gas flux measured by DOAS

Fig. 9
Fig. 9 Comparison between a HDR and other observations, including SO 2 flux, daily amplitude of continuous tremor (CT), lengths of the C-H and the H-S GNSS baselines, and daily sum of RMS of the long period tremor (LPT) amplitudes (product of the daily number of LPT and mean RMS per event) at Sunasenri station (Fig. 1a).The CT amplitude and LPT RMS are smoothed by seven-day median filter (black lines).The GNSS baseline lengths are smoothed by 2-month median filter (solid lines).identify the initiation timing of extension of the C-H baseline, contraction trend estimated during May 2020 to August 2021 was removed (red circles) from the daily baseline length (black circles).The red dashed line corresponds to the end of the 2019-2020 eruption.The black line corresponds to the timing of the crater depression.The solid red lines correspond to the 2021 eruptions.P1 to P4 indicate Phases 1 to 4. The blue highlighted zone corresponds to the May 2021 unrest.b The C-H baseline length around the timing of the baseline extension (May 2021 to December 2021)

Fig. 10
Fig. 10 Comparison between vent degassing and VLP activity.The daily RMS of the VLP is compared with a the H 2 O flux estimated from the thermal observation and b SO 2 flux measured by DOAS.The black, red, and blue symbols indicate values in Phases 2, 3, and 4, respectively