Continuous monitoring of the 2015–2018 Nevado del Ruiz activity, Colombia, using satellite infrared images and local infrasound records

Nevado del Ruiz Volcano (NRV) had a phreatomagmatic eruption in 1985. The eruption partially melted the volcano’s ice cap leading to floods and lahars flowing down to nearby towns, which killed at least 25,000 people. This event has raised particular importance of monitoring activity including small eruptions at ice-capped high-altitude volcanoes. However, the high altitude makes it difficult to maintain monitoring stations near the summit crater. Moreover, the visibility of the summit area is frequently prevented by clouds. In this paper, we report the results of a feasibility study for detecting thermal anomalies and small eruptions using satellite thermal remote sensing and ground-based infrasound technique. We newly included South and Central America to the observation areas of the near-real-time monitoring system of the active volcanoes, which uses infrared images from satellites. We also operated three infrasound stations in the distances of 4–6 km from the active crater. Each of the stations consisted of a pair of infrasound sensors, and a cross-correlation technique was applied. The thermal and infrasound data acquisition started in August 2015 and December 2016, respectively, and recorded the recent dome-forming activity of NRV. We proposed parameters representing the visibility of the thermal anomalies and infrasound signals. These parameters are useful for monitoring because the severe weather condition at NRV frequently prevents signal detections. We discussed the detected thermal anomalies and infrasound signals in comparison with their visibilities and the changes in the volcanic activity of NRV reported by the local observatory. The thermal anomaly and infrasound detections were consistent with the high eruptive activity occurring at the NRV from October 2015 to May 2017 and its subsequent decline. Within the active period, there were breaks in the detections of thermal anomaly and infrasound. The visibility analyses allowed us to interpret the breaks as a result of bad weather conditions and to distinguish them from the confirmed low-activity periods after May 2017.


Introduction
Monitoring of high-altitude volcanoes in tropical zones is very difficult, which is a serious issue in the mitigation of volcanic disasters in these areas. Thermal and infrasound monitoring techniques are widely used tools for tracking and quantifying the surface activity of volcanoes (e.g., Kaneko et al. 2010aKaneko et al. , b, 2019Ripepe et al. 2018;Laiolo et al. 2019).
Thermal anomalies are features that provide information about the surface thermal state of the volcanoes. Anomalous changes in the surface temperature are associated with variations of heat flux from the underground Open Access *Correspondence: lmcastano@sgc.gov.co 1 Servicio Geológico Colombiano, Bogota, Colombia Full list of author information is available at the end of the article to the surface (Oppenheimer 1998). Thermal variations at volcanoes can indicate unrest, eruptive crisis, ongoing eruptive processes, or precursors to more energetic explosive/effusive activity (Dean and Dehn 2015). Thus, the continuous recording of thermal data for different eruptive types, stages, or cycles in the volcano is very useful for distinguishing clearly the cause of the anomaly.
Satellite-based infrared sensors can detect volcanic thermal features. These instruments are useful tools for monitoring volcanoes that are active or potentially active but are poorly instrumented (e.g., Dean et al. 1998;Kaneko et al. 2010a). Because satellite sensors make daily or semi-daily observations of a region with moderate to high resolutions, a number of near-real-time systems based on automated algorithms have been implemented (Harris et al. 2001;Dehn et al. 2002;Wright et al. 2002Wright et al. , 2004Coppola et al. 2009Coppola et al. , 2016. REALVOLC (Kaneko et al. 2002(Kaneko et al. , 2010a is one of the automated monitoring systems for tracking thermal anomaly using satellite data from Moderate Resolution Imaging Spectroradiometer (MODIS) and Multi-functional Transport Satellite (MTSAT) sensors. This system has been in operation since 2010, monitoring 147 active volcanoes in East Asia, and it has captured the eruptive sequences of Mt. Asama (Japan) and Sarychev Peak (Russia) volcanoes (Kaneko et al. 2006(Kaneko et al. , 2010a. Infrasound signals from large eruptions are detected worldwide and are used for global monitoring of eruptions, including those at remote volcanoes. The pioneer works are Passechnik (1959) at the Bezymyanny volcano (Kamchatka, Russia) and Goerke et al. (1965) at the Agung volcano (Bali, Indonesia). In recent years, highperformance infrasound stations are operated as a part of the International Monitoring System (IMS) for verification of the Comprehensive Nuclear-Test-Ban Treaty (CTBT) and are also used for monitoring and studying volcanic eruptions (Evers and Haak 2005;Le Pichon et al. 2010;Matoza et al. 2017;Marchetti et al. 2019).
Local infrasound observation is necessary for monitoring small eruptions and gas emissions (Johnson et al. 2003;Lopez et al. 2013;Ulivieri et al. 2013), and for quick detection of eruptions (Ripepe et al. 2018). The wind noise problem is particularly severe for the detection of signals from small eruptions at high altitudes. Recording with multiple sensors in an array effectively distinguishes signals from noise (Ripepe and Marchetti 2002;Matoza et al. 2007). Most of the array observations were made temporally, while permanent infrasonic arrays are operated only at a limited number of well-monitored volcanos such as Stromboli and Etna in Italy (Ripepe et al. 2007;Ulivieri et al. 2013) and Kilauea in Hawaii (Thelen and Cooper 2014). An alternative method to distinguish infrasonic signals from wind noise was to make cross-correlation analyses between signals from an infrasound sensor and a collocated seismometer (Ichihara et al. 2012). The method was applied at many volcanoes with some improved analytical techniques (e.g., Cannata et al. 2013;McKee et al. 2018). However, the method does not work when a volcano is seismically very active (Ichihara et al. 2012).
In this study, we apply thermal and infrasonic methods designed for continuous monitoring for the 2015-2018 activity of NRV. A phreatomagmatic eruption of NRV in 1985 caused one of the most significant volcanic disasters in history. The eruption partially melted the volcano's ice cap leading to floods and lahars flowing down to nearby towns, which killed at least 25,000 people. This event has raised the particular importance of monitoring ice-capped high-altitude volcanoes. However, the high altitude makes it difficult to maintain stations near the summit crater. Moreover, the visibility of the summit area is frequently prevented by clouds. The satellite thermal monitoring system is operated since August 2015, and the new infrasonic stations are since December 2016. They recorded data during the recent dome-forming activity and a series of small eruptions. Because the activity was also routinely reported by the local observatory of Colombian Geological Survey (SGC-OVSM) with multiparametric monitoring data, we had a rare opportunity to make a feasibility study of detecting the surface activity by the thermal and infrasonic methods at a high-altitude volcano in a tropical region.

Recent activity of Nevado del Ruiz Volcano (NRV)
NRV is one of the most active volcanoes in Colombia with a height of 5321 m above sea level locating in the Cordillera Central de Colombia (4° 53′ 43″ N; 75° 19′ 21″ W), approximately 140 km North-West of Bogotá and 28 km South-East of Manizales city (Fig. 1a). It is an andesitic volcano with recent activity at least in the last ≈10 kyrs mainly explosive in nature (Martínez et al. 2014;Ceballos-Hernández et al. 2019 (Londono 2016), with surface manifestations at the end of 2010. Temporal changes in seismicity, deformation, and geochemistry (SO 2 release mainly) have been observed since then, with changes in the dynamics of the current eruptive episode (Castaño-López et al. 2017). Two small phreatomagmatic eruptions took place in May and June 2012 (Martínez et al. 2012). Later on, in October 2015, an effusive eruption occurred, with the emplacement of a small lava dome at the bottom of the active crater. Small Vulcanian eruptions have been recorded frequently during these eruptive periods (Londono and Galvis 2018). Figure 2 summarizes some of the parameters that SGC used for monitoring the activity of NRV (SGC Servicio Geológico 2015, 2016. Londono and Galvis (2018) calculated the volume of ash emission for selected events by 1D modeling of a volcanic plume using the observed column height from ground photographic images and the eruption duration obtained from the seismic signal synchronized with the images. The evaluated volumes are shown in Fig. 2c with daily counts of ash emission events. The possible volume of the extruded lava dome (Fig. 2e) was calculated from 57 available satellite images from TerraSAR-X (TSX)/TanDEM-X (TDX) radar satellites (see Appendix A). Combining Fig. 2 and reports from the volcano monitoring (e.g., SGC Servicio Geológico 2015, 2016, we have defined three stages of the volcanic activity for the period from November 2014 to December 2018. Stage I, from November 2014 to October 2015, corresponds to a period of magma ascent with high seismicity (Fig. 2a, b) and volcanic tremor (continuous and pulsating), edifice deformation, ash emissions (Fig. 2c), and SO 2 release ( Fig. 2d) (Lundgren et al. 2015;SGC Servicio Geológico 2015;Londono 2016). At the end of Stage I, a lava dome started emplacement at the bottom of the crater (Fig. 2e). Stage II, from October The cumulative numbers of volcano-tectonic (VT) earthquakes (red), long-period (LP) or very-long-period (VLP) events (blue), and hybrid (HB) ones (green). b The cumulative energy of VT (red) and HB (green) events. c Ash volume for the reported eruptions (triangles) obtained by (Londono and Galvis, 2018) and daily counts of ash emission events (bars) from SGC-OVSM catalog. d The daily maximum SO 2 flux (blue bars) and its cumulative values (orange) 2015 to May 2017, corresponds to a dome growth period with a new possible input of magma to the initial volume ( Fig. 2e), accompanying the significant increase of ash emissions (Fig. 2c), seismicity (Fig. 2a, b), SO 2 (Fig. 2d), and volcanic tremor episodes with moderate deformation (SGC Servicio Geológico 2015, 2016. Stage III, from May 2017 to December 2018, corresponds to final emplacement and subsequent lava dome evolution, which accompanied minor ash emissions (Fig. 2c), continuous and variable SO 2 release (Fig. 2d), low seismicity levels (Fig. 2a, b) and very low deformation. We will compare such division of Stages I-III with the detected thermal anomalies and infrasound data.

Thermal anomalies
REALVOLC is a near-real-time monitoring system used in active volcanoes to track the volcanic surface thermal state and evaluate their level of activity. Since August 2015, the REALVOLC system has started covering the Central and South America regions, with more than 80 active volcanoes in the Andes mountain range monitored. In Colombia, 11 out of the 23 active volcanoes have been included in the web-based system (http://vrsse rv2.eri.u-tokyo .ac.jp/index .html).
For NVR and other volcanoes, the system uses level 1b nighttime MODIS infrared images. These images are acquired by sensors mounted on the NASA Terra and Aqua satellites with a horizontal spatial resolution of 1 km per pixel. MODIS data are downloaded from the NASA LAADS web service and each image or scene is trimmed at 101 × 101 pixels as a preprocessing area. The analyzed area is 6 × 6 km centered on the volcano summit ( Fig. 3a). We process and analyze, through an automatic main routine, the calibrated radiances from the nighttime MODIS 21 and 31 bands with the central wavelength at 3.9 μm and 11 μm, respectively (Kaneko et al. 2006(Kaneko et al. , 2010a. Calibrated radiances at the MODIS 3.9 μm and 11 μm bands were converted to brightness temperatures using the Planck function (Rothery et al. 1988). The band 21 is used for detecting thermal anomalies on volcanoes, while band 31 is for inferring the background temperature of the region (Wooster and Rothery 1997). The obtained data from both bands are plotted together to show the thermal anomaly compared with the background temperature (Fig. 3b).
In the analysis of background temperature level, we used time-series data from Nevado del Tolima volcano (NTV) approximately 26 km at the South-Southwest of NRV ( Fig. 1). This analysis was done to establish a baseline related to the regional climatology, as NTV (height 5220 m a.s.l) has topographic and glacial characteristics similar to NRV, while it showed low volcanic activity levels during the studied period.
Thermal radiance from the NRV is very low when the target area is covered by clouds. To eliminate these unreliable low-temperature values due to the poor visibility of the ground surface, a statistical data filtering method was applied. For a sliding time window of ± 10 days centered at each day, we selected the maximum 10% of Band 31 values and took their mean. This method gave us the baseline of the temperature, excluding the data affected by clouds. The threshold to accept a value as a data point was that both values of Bands 31 and 21 were larger than the baseline value. The fraction of the accepted data points in sliding 4-day windows was used as the visibility of the thermal anomaly of the period, which we refer to ThVis. The length of 4 days was used because the data intervals are approximately 1 day (one or two points during the nighttime every day) and the cloud condition can change in less than a week. It is noted that ThVis thoroughly represents the weather condition and is independent of the volcanic activity.

Infrasound
Some infrasound stations had been operated at NRV, each of which had a single infrasound sensor collocated with a seismometer. However, we could not resolve infrasound signals covered with wind noise. The crosscorrelation method using an infrasound sensor and a collocated seismometer (Ichihara et al. 2012) was not useful at NRV due to the vigorous seismic activity. Therefore, we have newly installed three infrasound stations (ACOLM, ACINM, and ACRUM stations, Fig. 1a), each of which has a pair of infrasound sensors separated by 5-8 m.
The sensor pair consists of a Chaparral 60Vx (frequency range 0.1-200 Hz, full scale 720 Pa, noise level 2 mPa RMS in 0.5-2 Hz) and a Hakusan SI104 (frequency range 0.3-2000 Hz, full scale 2 kPa, noise level ~ 0.1 peak to peak Pa). Although the latter sensor has a larger instrumental noise, it has a sensitivity below 0.1 Pa. A Guralp DM26 digitizer records them at 100 Hz in miniSEED format. Continuous records were available from December 2016. We mainly used ACOLM station because it is the closest to the active crater (4.2 km) and had the best signal quality. The analysis of infrasound activity was done from December 2016 to November 2018.
We performed the cross-correlation analysis for the two collocated infrasound sensors ('C' and 'H' in Fig. 1a) to distinguish infrasound signals from wind noise. The detail of the method is explained in Appendix B. The concept of the method is similar to that of Ichihara et al. (2012), who analyzed the records of an infrasound sensor and a collocated seismometer. A frequency band from 1.5 to 8 Hz was used to guarantee that the distance between Fig. 3 An example of automatized analyses by REALVOLC. a The target area at NRV for the thermal analyses. b Thermal radiance as a function of Julian day. The radiance of a pixel in the 3.9-μm band (red crosses) responds to hotspots in the pixel, while that in the 11-μm band (green crosses) represents the mean temperature. We use the 3.9-μm band to detect thermal anomaly and the 11-μm band to estimate the background temperature the sensors was greater than the correlation length of wind noise (Shields, 2005) and much shorter than the wavelengths of infrasound signals. The cross-correlation coefficient of the filtered signals from the sensor pair was calculated using a 5-s sliding window centered at time t , with an overlap of 4 s. It is denoted as CC(X C , X H ; τ ; t) in Eq. (B.3), where X C and X H are the filtered data from the sensors 'C' and 'H' , respectively, and τ is the time lag of X H to X C . A daily plot of the cross-correlation coefficient as functions of the t and τ was analyzed by performing a manual search in every 1-h window to detect cross-correlation patterns associated with infrasound (See Additional file 1).
We considered the average over 20 s of t . A theoretical time lag τ of an infrasound arrival to the sensor pairs from the crater is approximately zero because the two sensors are close and aligned perpendicular to the direction of the crater. However, the direction of wave propagation can fluctuate due to the wind. Therefore, we take the maximum in −0.02 < τ < 0.02. The resultant value of CC is represented as CC(X C , X H ; 0; t).
Like the thermal and visual observations are prevented by clouds, the infrasound signal detection is influenced by wind noise. We define the 'visibility' of infrasound signals as InfVis. To evaluate InfVis, we calculated the mean-square values, MS(X n ; t), of the data from each sensor (n is either C or H) in the same frequency band (1.5-8 Hz) and time windows (20 s) centered at t . It includes powers of infrasound signal ( A n ), wind noise ( W n ), and instrumental noise ( N n ), where the subscript indicates the sensor 'C' or 'H' . If MS becomes large due to an increase of power of A n , CC also increases. When only MS increases with CC remaining small, on the contrary, W n or N n is so large that A n is obscured if any. See Appendix B for a further explanation. We compared the maximum value of CC(X C , X H ; 0; t) (MaxCrr) and the median value MS(X C ; t) (MedMS) in every 200-s time window. The time lengths of 20 s and 200 s were chosen, referring to the typical lengths of observed infrasound signals at eruptions. As MedMS mainly represents the noise level, the daily fraction of time windows in which MedMS is smaller than a threshold is defined as InfVis. The threshold value is determined by the following analysis.
Detectability analyses for infrasound signals were performed for all three stations by comparing the infrasound correlation patterns with the eruption event catalog. We made the catalog combining the monitoring data (SGC Servicio Geológico 2016, including seismic and photographic records as well as eyewitness reports (See Additional file 2: Table S2.1). It includes 311 confirmed eruptions, of which ash plumes are visually observed in 294 cases. The result of infrasound detection was also listed in the catalog (See Additional file 2: Tables S2.1 and S2.2). When we recognized the infrasound signal in the filtered wave trace as well as in the correlation pattern, we list the peak pressure value. The infrasound signal detection or no-detection was then compared with MedMS that represented the noise level.

Thermal anomalies
Figure 4a-c shows the results of the analyses for the thermal anomaly and its visibility (ThVis explained in the "Method" section) in the studied period (August 2015-November 2018). We applied the same method to the data of the reference volcano, NTV, and confirmed that the value in the Band 21 was rarely larger by 5 degrees than the baseline (Fig. 4c). Therefore, we regarded as a thermal anomaly when Band 21 showed a temperature larger by 5 degrees than the baseline. Out of the 8452 time windows at NRV, we accepted 1949 time windows (23%) as data points (the Band 21 values are larger than the baseline), 59 of which showed significant thermal anomaly.
NRV showed three periods of thermal anomalies with brightness temperature values ranging between 279 K and 320 K (Fig. 4b). Period 1 (P1) from September 2015 to June 2016 had the highest temperature (319.7 K) detected on December 31, 2015. In Period 2 (P2) from October 2016 to May 2017, we observed the maximum value of 315.5 K on March 4, 2017. The maximum value in Period 3 (P3) from September 2017 to March 2018 was 297.5 K on January 9. These periods of thermal anomalies were not observed at NTV (Fig. 4c).
The breaks between P1 and P2 and between P2 and P3 seem to correlate with the low thermal visibility (ThVis) at NRV (Fig. 4a). Therefore, it is not obvious whether the thermal activity was low in these breaks only from thermal data. We discuss this point later in "Discussion" section comparing the results with other monitoring parameters. Figure 5 shows two example events with infrasound signals recorded at ACOLM during the studied period. Although the signals are covered by wind noise in the raw data traces and are neither distinct in the filtered waveforms, the cross-correlation analyses between the collocated infrasound sensor pairs help to distinguish the infrasound from wind noise. Clear infrasound correlation patterns like those in Fig. 5 were observed in two periods (Table S2.1 in Additional File 2). The first was from December 2016 to May 2017, with acoustic signals associated with ash emissions lasting from 30 to 400 s, with peak pressures ranging from 0.06 to 2.5 Pa in the filtered waveforms. The second was from July to October 2018, with a maximum pressure value of 2.2 Pa.  Table S2.1) are marked, distinguishing the cases that the infrasound signals are noticeable as the cross-correlation pattern (detected) with magenta triangles and not (failed) with black crosses. The green triangles indicate the cases that have clear infrasonic cross-correlation patterns with the volcano's seismic activity but without the confirmed eruptions by cameras or external reports. (See Additional file 2: Table S2.2.) When MaxCrr is large that is the infrasound signal is dominant in the data, most of MedMS values are smaller than 10 −2 Pa 2 . It indicates that the average power of the infrasound signal in the time window is smaller than 10 −2 Pa 2 . When MedMS > 10 −2 Pa 2 , it is rare to have MaxCrr larger than 0.4 because wind noise dominates the infrasound signal in the time window. At ACINM (Fig. 6b) and ACRUM (Fig. 6c), the eruption signals are marked as 'detected' even with small values of MaxCrr. These signals were not clear but have been recognized helped by the detection at ACOLM.

Infrasound
The right column of Fig. 6d-f shows the numbers of the eruption events with infrasound detection (magenta) and failure (black) versus the MedMS in the time windows corresponding to the eruption events in the catalog. For the failure cases, we counted only the events that accompanied ash plumes higher than 200 m. sThe detection is very poor when the MedMS is higher than 10 −2 Pa 2 at all the stations. The detections are much less at ACINM and ACRUM than at ACOLM even when MedMS is small enough. We consider the reasons are (1) ACOLM is located closer to the crater than ACINM and ACRUM and (2) it is in the main downwind direction.
We evaluate the infrasound visibility (InfVis) by the daily fraction of the 200-s time windows that had MedMS < 10 −2 Pa 2 at each station ( Fig. 7b-d), similarly to the thermal visibility (ThVis) in Fig. 3. We also had wind data at the weather station close to ACOLM (Fig. 1a). We used the average wind speed data recorded with intervals about 10 min and calculated the daily fraction of the data smaller threshold values of 2, 4, and 8 m/s (Fig. 7a). The threshold value of 4 m/s gives a similar pattern to InfVis at ACOLM (Fig. 7b), confirming that InfVIs represents the wind noise condition. From December 2016 to May 2017, the InfVis was relatively high, and many signals are detected with eruption events. In both years of 2017 and 2018, the InfVis was low from June to September. In these periods, thermal visibility was also low (Fig. 4). We infer that few events were recognized because the seasonal weather condition caused poor visibility. On the other hand, from October 2017 to June 2018 and from October 2018 to November 2018 (marked by the blue arrows in Fig. 7), the InfVis was high at all the stations.

Comparison with the other monitoring parameters
The detected thermal anomalies and infrasound are compared with other parameters of the volcanic activity of NRV (Figs. 2 and 8a, b) with consideration of the seasonal change in their visibility.
The thermal analyses started at the end of Stage I when the lava dome emplacement began while the infrasound analyses started at the end Stage II during the lava dome growth (Fig. 2e).
The most significant thermal anomalies (P1 and P2) and infrasound signals are observed in Stage II. The maximum temperature during the studied period was recorded in P1 and the maximum peak pressure concurs with the beginning of P2. The high temperature may be associated with the extrusion of the lava dome. Such high increases in temperature have been observed at different volcanoes during lava dome extrusions (Wooster and Rothery 1997;Smith et al. 2011;Reaht et al. 2016). From the end of August 2016, the growth rate of the lava dome increased (Fig. 2e). It may have generated the second peak of the thermal anomalies (P2) and high level of infrasound activity.
From May 2017, the growth rate of the lava dome decreased, and the volume stabilized around 1 × 10 6 m 3 during Stage III (Fig. 2e). The P3 period of thermal anomalies was still observed in the first part of Stage III, but the temperature was lower than in P1 and P2. The P3 coincides with the confirmed low level of infrasound activity (the blue arrows in Figs. 7d and 8a). There were some registered eruptive events in this period but the number of clear ash emission events was only two (November 24 in 2017 and February 6 in 2018 as in Additional file 2: Table S2.1). We infer that when lava dome growth was stopping (Fig. 2e), the volcanic system remained open to allow gas or gas-andash releasing to the atmosphere without significant infrasound signals, and hot enough to produce thermal anomalies.
It is still uncertain whether the breaks in the thermal anomalies between P1, P2, and P3 were real or not. The visibility of infrasound was also low in the breaks due Fig. 7 The visibility of infrasound signals (InfVis) at NRV evaluated by the wind data at the OLLETA weather station close to ACOLM (a) and by the infrasound data at ACOLM (b), ACINM (c), and ACRUM (d). a The daily fraction of wind speed data smaller than the threshold values of 2 m/s (green), 4 m/s (blue), and 8 m/s (pink). b-d The InfVis is shown by blue, and the periods of missing data are shaded. The events to which the associated infrasound signals are detected and failed at the individual station are indicated with the magenta triangles and the black crosses, respectively (see Additional file 2: Table S2.1). The green circles are the unconfirmed events (see Additional file 2: Table S2.2). The blue arrows at the bottom show the periods of low infrasound activity, confirmed by the high InfVis to strong wind (Fig. 7). Therefore, the breaks may be due to the poor visibility (Fig. 4a).
During the second part of Stage III, no remarkable thermal anomalies were detected even with the sufficiently high visibility values (Fig. 4). A temporal increase in infrasound emissions was observed (Fig. 8a) at the end of Stage III (July to September). The lava dome volume did not increase but varied (Fig. 2e). Such variations may be associated with a new magma batch that ascended to shallower regions, so the intrusion did not affect considerably the temperature of the previous magma batch that formed the first lava dome. The permeability of the conduits may have changed with this new magma input (Pallister et al. 2013;Gaunt et al. 2016), sealing microfractures that partially pressurized the conduits, and allowing the generation of small explosions that exhibited infrasound signals. Figure 8c, d shows the number of icequakes (glacier earthquakes) recorded in the seismic stations BISM, RECM and ALFM (Fig. 1a) and air temperature time series of the OLLETA weather station (Fig. 1a). Icequakes are coseismic brittle fracture events within the Fig. 8 Comparison of infrasound and thermal data with surface-activity parameters recorded between September 1, 2014, and December 31, 2018. a Maximum pick pressure from the Chaparral sensor data at ACOLM station since December 2016. The periods shown by the blue arrows are the same as those in Fig. 7. b The same thermal anomaly data as in Fig. 4b. c Daily counts of icequakes at NRV (gray) and NTV (light blue). d Air temperature recorded by OLLETA weather station ice generated in glacier-covered areas (Ekström et al. 2003;Podolskiy and Walter 2016), which sometimes produce surface break-off events with infrasound signals (e.g., Preiswerk et al. 2016). There are multiple proposed mechanisms that cause icequakes especially on volcanoes covered by icecap, which include ice falling, ice surface cracking, resonant water-filled ice cavities, sudden changes in water flow rate (St. Lawrence and Qamar 1979;Métaxian 2003), stick-slip sliding (Allstadt and Malone 2014), and crevasse (Neave and Savage 1970). On the other hand, icequakes on volcanoes can be triggered by volcanic activity (Delgado et al. 2015) and/ or thermal stress (Lombardi et al. 2019). The icequakes induced by thermal stress happen when the air temperature decreases during clear nights (no clouds). The number of icequakes increases (Carmichael et al. 2012;Podolskiy et al. 2018;Zhang et al. 2019) possibly due to thermal contraction (caused by significant temperature fluctuation), which induces ice fractures (Lombardi et al. 2019). In Fig. 8b-d, it is possible to observe the increase in daily averaged air temperature and icequakes in synchronicity with thermal anomalies for P1, P2, and P3. Because the visibility of thermal anomalies (ThVis in Fig. 4) is better in clear nights, seasonal atmospheric variations may influence the occurrence of icequakes as well as the ThVis. Nevertheless, we cannot exclude the possibility that the recent surficial activity of NRV affects both of the thermal anomalies and icequakes occurrences. We can see that the number of icequakes suddenly increased in December 2014, and the peak of the icequake numbers in 2018 when the volcanic activity has declined is smaller than the previous 4 years. Also, the number of icequakes in these years at NRV is significantly larger than that at the less active neighbor volcano, NTV (light-blue color in Fig. 8c). A detailed work about this topic is needed to clarify the possible link between icequake and volcanic activity.

Infrasound monitoring for small eruptions at a high-altitude volcano
In this study, we installed a pair of infrasound sensors separated by less than 10 m at each station. The crosscorrelation analyses of the data from the sensor pair were useful for detecting weak infrasound signals, which a single sensor could not have resolved. On the other hand, while 101 confirmed eruptions were recorded by either of the infrasound stations, 32 events were detected simultaneously at all the stations (see Additional file 2: Table S2.1). The ACOLM station detected more events than the others, of which reason we consider that ACOLM was the closest to the active vent and located in the main downwind direction as mentioned in the previous section. Although the atmospheric propagation effect is considered significant at larger distances, Lacanna et al. (2014) show that propagating infrasonic power is directional as a result of combined effect of wind and topography even within kilometers. There were also cases that some stations missed the signals due to locally higher noise levels.
The detection would be better if a larger number of infrasound sensors are installed in an array at a station. However, when we have a limited number of sensors, we consider it challenging to concentrate all the sensors at a single station. We propose that the current method, which consists of sensor pairs at multiple stations, is useful for monitoring a high-altitude volcano. The sensor distance is small enough to connect them to a single recorder and a power source so that the costs for installation and maintenance are saved. On the other hand, the current system fails to resolve the infrasound source directions. The main disadvantage is that the sensor pairs are aligned perpendicular to the direction of the active vent. It is better to install them in the line of the vent direction. Even if the sensor separation is as small as 10 m, an infrasound wave propagating from the vent has a time delay larger than 0.02 s, which we could distinguish by the 100-Hz recording.

Concluding remarks
We made a feasibility study of detecting small eruptions and thermal anomalies of NRV. We analyzed thermal data from satellite remote sensors and infrasound data from newly installed three stations for 2 years and compared the results with the volcano monitoring data by the local observatory (SGC-OVSM). Detection of infrasound signals from small eruptions has been significantly improved by installing a pair of infrasound sensors separated by 5-8 m. Of the three infrasound stations, the closest one located 4.2 km from the crater to the east made significantly better detection than the others.
Spite of thermal anomalies and infrasound signals were not always detectable over all the studied period, due mainly to weather conditions, we could detect a good number of them and quantify their visibility. Based on the combined analysis using the satellite infrared and infrasound observations, we infer that the activity, which started in November 2014, entered into the main active phases between October 2015 and May 2017, and lowered later on. Within the active period, there were gaps in the detections of thermal anomaly and infrasound. The visibility analyses allowed us to interpret the breaks as atmospheric effects and to distinguish them from the confirmed low-activity periods after May 2017.
Satellite observation of thermal anomalies and infrasound technique is required for monitoring the surface activity of high-altitude volcanoes in tropical zones, which frequently have clouds and strong wind. The visibility analyses are essential for interpreting the data. The methods and results shown in this study will be useful for future works in this field.

Supplementary information
Supplementary information accompanies this paper at https ://doi. org/10.1186/s4062 3-020-01197 -z.   The data recorded by a sensor n(n = CorH) is represented by X n . We represent the cross-correlation function between X C and X H by R(X C , X H ; τ ) , where τ is the time delay of X H to X C . The cross-correlation coefficient, CC(X C , X H ; τ ) is where P(X n ) = R(X n , X n ; 0) , which is the power of X n . We assume the data X n consists of pressure fluctuation from infrasound signals, A n , seismic-to-infrasonic converted wave, S n , wind noise, W n , and instrumental noise, N n . Namely, The component S n is usually negligible unless the ground velocity is very large. We have confirmed using some earthquake events that the infrasound stations used in this study recorded pressure oscillation less than 1 Pa for the vertical ground velocity of around 1 mm/s, which is in the typical range of the local atmospheric response to the ground motion (Ichihara et al. 2012;Watada et al., 2006). Because the ground velocity of seismic waves related to eruptions at NRV is too small to be recorded by the infrasonic stations, we omit S n in the following discussion.

Abbreviations
In principle, there is no correlation between each other of A n , W n , and N n . It is also confirmed that the correlation of the instrumental noise of the collocated two sensors is negligible. Then, the cross-correlation coefficient between the sensor 'C' and 'H' is approximated by Shields (2005) investigated the correlation length of wind noise using a microphone array and showed that the correlation decays with distance, x , as where f is the frequency and v is the wind speed. Substituting (B.4) to (B.3) and assuming P(W C ) ∼ P(W H ) , we obtain In order to detect the infrasound signal, the second term of the numerator must be small enough compared with the first term. The second term is reduced if 3.2xf v > 1 . Considering x > 5 between the two sensors and assuming v < 20 for most of the time, we used the frequency range of f > 1.5 Hz.
(1) CC(X C , X H ; τ ) = R(X C ,X H ;τ ) √ P(X C ) √ P(X H ) , (2) X n = A n + S n + W n + N n (4) CC(X C , X H ; τ ) = R(W C ,W H ;τ ) CC(X C , X H ; τ ) ∼ The time lag τ which makes R(A C , A H ; τ ) the largest is expected to be close to zero, because the sensor pairs at the individual stations are close and aligned perpendicular to the direction of the crater (Fig. 1). However, the propagation direction can be locally fluctuate depending on the wind conditions, we evaluate CC(X C , X H ; τ ) by the mean value for −0.02 < τ < 0.02 s. Because the CC changes sharply with τ when the wavelength of the infrasound is not long enough compared with the sensor separation distance, we limit the frequency used in the analyses below 8 Hz, which corresponds to the wavelength of 40 m for the sound speed of 320 m/s. Equation (5) also indicates that when P(W C ) ≫ R(A C , A H ; 0) ∼ P(A C ) , CC is significantly reduced by the wind noise. In other words, we may detect an infrasonic signal when it arrives in time windows with small P(W C ) . Therefore, we define the 'acoustic visibility' by the probability of the time window in which P(W C ) is smaller than a threshold. The threshold depends on the strength of the infrasound signals of interest.