Diffuse CO₂ degassing precursors of the January 2020 eruption of Taal volcano, Philippines

Abstract

On January 12, 2020, Taal volcano in Philippines erupted, 43 years after its previous eruption in 1977. This eruption was preceded by diffuse CO₂ degassing precursory signals. Significant temporal variations in diffuse CO₂ emission from Taal Main Crater Lake (TMLC) were observed across the ~ 12 years reaching high CO₂ degassing rates in 2011 and 2017, with values typical of plume degassing volcanoes. In addition to these CO₂ surveys at the TCML, soil CO₂ efflux continuous monitoring was implemented at Taal volcano since 2016 and a clear increasing trend of the soil CO₂ efflux in 2017 was observed. These geochemical observations are most simply explained by magma recharge to the system, and represent the earliest warning precursor signals to the January 2020 eruptive activity.

Introduction

Taal Volcano Island is the scenario of powerful eruptions and is the largest volcanic threat to the Philippines. The thirty-three recorded eruptions of Taal Volcano Island from 1572 to 1977 include phreatic to phreatomagmatic eruptions. Six of the 33 known eruptions since 1572 have resulted in fatalities^1,2 and today several million people live within a 20-km radius. Because of these facts, Taal Volcano Island (Fig. 1) was one of the 16 Decade Volcanoes³ identified by the International Association of Volcanology and Chemistry of the Earth's Interior (IAVCEI) as worthy of particular study in light of their history of large, destructive eruptions and proximity to populated areas after the United Nations General Assembly designated the 1990s as the International Decade for Natural Disaster Reduction (IDNDR).

Figure 1

(a) Geographical location of Taal Volcano Island, Philippines. (b) Aerial view of the volcano taken from the NE showing the TMCL. (c) Shaded relief map of Taal Volcano Island showing the Taal Main Crater Lake (TMCL) and the location of the automatic continuous geochemical monitoring station (blue star). The map was constructed with the software Surfer version 8.00 surface mapping system (Golden Software, Inc; https://www.goldensoftware.com/products/surfer).

Following this international awareness, a collaborative geochemical monitoring research program between Philippine and Spanish scientists was established to contribute to the strengthening of volcanic surveillance of Philippine volcanoes. This collaborative research was focused mainly in the monitoring of diffuse CO₂ degassing since it is the CO₂ major gas component, beside water vapor, in both volcanic-hydrothermal fluids and magmas. It is also a good tracer deep of sub-surface magma degassing, since its low solubility in silicate melts at low to moderate pressure favors its early exsolution^4,5,6. Volcanogenic CO₂ is released not only through preferential degassing routes in volcanoes such as fumaroles and plumes, but it could partially also percolate through the entire volcanic edifice and released to the atmosphere in a diffuse or “silent” mode”^7,8,9. In addition, large quantities of thermal energy are released by volcanoes through its diffuse CO₂ emission^10,11,12 and this type of degassing can be also a significant contributor to the subaerial global volcanic CO₂ degassing^13,14,15,16.

Diffuse volcanic degassing disturbs the chemical and isotopic composition of the soil-air and water–air interfaces at the surface environment of the volcano, producing enrichments not only of CO₂ but also of He, H₂ and other tracer gases^17,18,19. One of the first studies of diffuse degassing on volcanoes was about continuous soil gas H₂ monitoring at Mount St. Helens¹⁷. During the last 25 years numerous gas geochemical studies have highlighted the importance of this type of degassing in volcanic systems^{7,8,9,20,21,22,23,24,25,26} and its great use to strengthen the geochemical monitoring program for volcanic surveillance^{27,28,29,30,31,32,33}, particularly at those volcanic areas where visible volcanic gas emissions (plume, fumaroles, etc.) either are scarce or do not exist^34,35,36,. However, the detection of diffuse CO₂ degassing anomalies prior to volcanic eruptions reported are very scarce^{34,37,38,39,40}. Hydrothermal fluids, a mixture of seawater, volcanic water and meteoric water⁴¹, feed the surface discharges of Taal volcano and produced strong hydrothermally altered areas exhibiting solfatara, fumaroles, hot springs, and gas bubbling in the TMCL. Taal volcano has suffered frequent periods of unrest since the eruption in 1977, characterized by increases in the seismic activity, ground deformation and gas emissions. Magmatic intrusions have been the most plausible mechanism to explain these unrests^30,42,43, although other authors have pointed out that magma intrusions were very unlikely after 1994⁴⁴. On 12 January 2020, a volcanic eruption occurred in the main crater of Taal volcano. The eruption was characterized first by a phreatic-phreatomagmatic style, producing a giant plume of volcanic ash up to ~ 15 km in the atmosphere⁴⁵, and ended with a less explosive eruption characterized by the occurrence of lava fountains. The eruption ejected juvenile products, with evidence of magma mingling⁴⁶ and represented a great impact to the Philippines, as around half a million people were directly affected by the event, producing the loss of ~ 69 M$ worth of damage to infrastructure and agriculture⁴⁷. The acidity of the TMLC waters (pH ~ 3) allows the emission of big amounts of CO₂ to the atmosphere, as at low pH values, the water of the lake reduces dramatically its ability to dissolve acid gas species as CO₂. Thus, monitoring CO₂ emission through the water surface is an important monitoring tool to detect early warning signals of future volcanic unrest episodes. We report here the earliest precursory signal of this volcanic eruption obtained after hundreds of diffuse CO₂ emission measurements covering the entire surface of the main crater lake and through a continuous monitoring of this gas in a single observation site.

Results

About 2630 CO₂ efflux measurements have been performed covering homogenously the 1.2 km² of the TMCL across the ~ 12 years through 19 surveys (an average of ~ 138 measurements per survey, what means ~ 115 measurements/km²) showing a wide range of values from > 0.5 g m⁻² d⁻¹ up to 84,902 g m⁻² d⁻¹. Table 1 summarizes the diffuse CO₂ emission rate (the amount of CO₂ that was being released to the atmosphere through the water surface of the TMCL at the time of the survey) observed in the period 2008–2018. Figure 2 shows the spatial distribution of the CO₂ efflux values in the period 2016–2018, where an important increase of the diffuse CO₂ emission values can be observed, reaching a relative maximum in November 2017. Statistical-graphical analysis of CO₂ efflux data of each survey at the TMCL has shown two geochemical populations (background and peak) suggesting the occurrence of either two different CO₂ sources or degassing dynamics, i.e. advection versus diffusion weight factor on the CO₂ degassing processes as it has been observed in other volcanic lakes⁴⁸. The background and peak populations of each survey are characterized by relatively low and high CO₂ efflux values and their average mean values are 643 g m⁻² d⁻¹ and 3707 g m⁻² d⁻¹, respectively.

Table 1 Summary of the diffuse CO₂ emission at TMCL in the period 2008–2018.

Figure 2

Spatial–temporal variations of CO₂ efflux measurements at the surface of the TMCL from 2016 to 2018.

The background diffuse CO₂ emission from TMCL has been estimated by multiplying the CO₂ efflux values of the background mean ($\overline{{\text{x}}}$) and range (± 1σ) population times the surveyed area. Therefore, the estimated background diffuse CO₂ degassing from TMCL shows an average of 782 t d⁻¹ with a ± 1σ range values of 1288 and 508 t d⁻¹. This background or baseline emission value, established after ~ 12 years of monitoring this geochemical parameter, is in good agreement with the background diffuse CO₂ emission observed during the period 2008–2010, with a range value between 506 ± 15 and 947 ± 22 t d⁻¹. The observed relatively high and anomalous diffuse CO₂ emission rate along the ~ 12 years reached values of 4670 ± 159 t d⁻¹ on March 24, 2011, and 3858 ± 584 t d⁻¹ on November 11, 2017, which were higher than $\overline{{\text{x}}}$ + 2σ values (2645 t d⁻¹) (Fig. 3a). It is worth noting that the $\overline{{\text{x}}}$ + 1σ value was exceeded in 15 November 2014 and 15 April 2015. Figure 3b depicts the correlation between the mean of the two diffuse CO₂ efflux geochemical populations (background and peak) and the CO₂ emission rate and reveals a higher influence of the peak population in the CO₂ emission rate with its increasing.

Figure 3

(a) Temporal variations of the diffuse CO₂ degassing rate at the TMCL from 2008 to 2018 (black solid circles and line) and number of volcanic earthquakes per week at Taal Volcano (grey columns) constructed following different published data: from 2008 to 2012²⁷ and from 2013 to 2020³⁹. (b) Correlation between the mean of the two diffuse CO₂ efflux geochemical populations (background and peak) and the CO₂ emission rate.

In order to provide an additional geochemical tool to improve the volcanic surveillance program of Taal, an automatic geochemical station able to measure soil CO₂ efflux with an hourly frequency, was installed at Daang Kastila (DAK), near the main fissures field on the northern flank of the volcano (Fig. 1). The complete time series measured by the geochemical station is composed of 22,414 valid measurements of soil CO₂ efflux, wind speed, air humidity and temperature, barometric pressure and soil temperature and moisture, from 25 January 2016 to 31 August 2019. Due to instrumental and power supply problems, time series has a 28.9% of missing data. A probability plot of the data allows us to distinguish two main log-normal geochemical populations: background (70% of the data) and peak (3.8% of the data), with 0.14 kg m⁻² d⁻¹ and 0.50 kg m⁻² d⁻¹ means values, respectively. The average value of the soil CO₂ efflux data showed oscillations around background values until 14 March 2017. Since that date at 22:00 h, the station measured a sharp increase of soil CO₂ efflux from ~ 0.1 up to 1.1 kg m⁻² d⁻¹ in 9 h and continued to show a sustained increase in time up to 2.9 kg m⁻² d⁻¹ in 2 November 2017, that represent the main long-term variation of the soil CO₂ emission time series (Fig. 4). Additionally, the diffuse CO₂ emission survey carried out on 16 March 2017 at TMCL, showed a relative peak value of 1763 t d⁻¹, and was followed by the second maximum value of the diffuse CO₂ emission survey series (3858 t d⁻¹) measured in 16 November 2017, a couple of weeks after the observed maximum value by the automatic geochemical station. One year later, a diffuse CO₂ emission survey carried out on 22 November 2018 at TMCL, showed a slightly lower value (3050 t d⁻¹). The emission rates measured after 2016 were ~ 3.7 times higher in average than the estimated background diffuse CO₂ emission from TMCL. Similar increases in the CO₂ released by the TMCL were reported in 2010–2011³⁰ and observed in 2014–2015.

Figure 4

One week moving average of the temporal variations of the soil CO₂ efflux values measured at the automatic geochemical station (red solid line) and diffuse CO₂ degassing rate measured at the TMCL (black solid circles and line), in the period January 2016–August 2019. Grey columns show the number of earthquakes per week³⁹. Dashed horizontal line depicts the log-normal background population ($\overline{{\text{x}}}$) and dotted horizontal lines depict the background the range ($\overline{{\text{x}}}$ ± 1σ) and ($\overline{{\text{x}}} \pm 2\sigma$) obtained after a probability plot analysis of the soil CO₂ efflux values measured at the automatic geochemical station.

Discussion

The biogenic contribution to the background population at the TMCL seems to be negligible, as the CO₂ efflux values measured in the Taal Caldera lake reported previously by other authors⁴⁹, where the observed background and peak mean population values (1.3 g m⁻² d⁻¹ and 10.5 g m⁻² d⁻¹, respectively), were two orders of magnitude lower than those observed at the TMCL. Deep-seated CO₂ degassing has existed from TMCL during the study period and the occurrence of two major populations for the diffuse CO₂ emission at the TMCL could be closely related to differences of the gas transport mechanism (diffusion and advection) which might be controlled by permeability variations in the volcanic system. The spatial distribution maps of diffuse CO₂ emission at the TMCL show significant spatial–temporal variations of CO₂ efflux measurements for the period 2008–2012³⁰, and such variations can also be observed for the last 5 surveys performed during the period 2016–2018 (Fig. 2). The observed diffuse CO₂ emission values from TMLC across the ~ 12 years, reaching high CO₂ degassing rates in 2011 and 2017, are typical of plume degassing volcanoes¹⁴ and highly active volcanic systems^14,33.

The maximum diffuse CO₂ emission rates across the ~ 12 years (24 March 2011 and 11 November 2017, Fig. 3a) are not only higher than the + 1σ background range values, but also higher than the + 2σ ones (~ 2518 t d⁻¹) suggesting higher injection rates of magmatic fluids into the volcanic-hydrothermal system of Taal. The 2010–2011 unrest phase was characterized by significant increases in CO₂ emission³⁰, but the maximum CO₂ emission rate measured in that period occurred 2 months before the strongest seismic activity recorded during the unrest period. Other episodes of magmatic fluids injection into the volcanic-hydrothermal system likely occurred between 2014 and 2015⁴³, and the main increase in seismic activity that generated this magma movement was again preceded by 2–3 months by the increase in CO₂ emission observed in November 2014. The largest time gap between the observed geochemical anomaly (CO₂ emission rate) and the seismic activity took place in the third magmatic intrusion that led to the 12 January 2020 eruption, as the maximum CO₂ emission rate released through the water surface of TMCL occurred in November 2017 or before (Fig. 3a), roughly 17 months before the significant increase recorded in the seismic activity. This increased difference in the time gap between the observed CO₂ emission rate and the start of anomalous seismic activity suggest a deeper source for the third magmatic intrusion that cause the last CO₂ emission rate anomaly due to decoupling between the amount of volatile exsolution and magma rise. After November 2017, the temporal variations of the soil CO₂ efflux values measured at the automatic geochemical station showed an important decrease until it recuperated background values (Fig. 4), likely due to the degassing of the magma that intruded deep at the beginning of 2017.

The inspection of Fig. 3b shows that peak values of the CO₂ effluxes control in greater proportion the CO₂ emission rate. Indeed, the slope of the correlation between the mean peak geochemical populations and the CO₂ emission rate is 5.39 times greater than the background slope. This observation highlight the importance of monitoring the amount of CO₂ emitted by the lake taking into account the whole surface and not only one observation site, as the selected observation site can be located in an area that exhibits background emission values during the phase of unrest. In fact, and as it is observed in Fig. 2 and has been observed by other authors³⁰, the location of CO₂ emission anomalies varies greatly in TMCL. The highest diffuse emission rates (> 2000 t d⁻¹) measured during the study period correspond to the surveys in which the highest peak values were measured (> 3000 g m⁻² d⁻¹), with background emission values < 2000 g m⁻² d⁻¹. This fact can only be explained by the injection of hot magmatic fluids from gas-rich magma into the hydrothermal system and subsequent escape towards the surface through the TMCL.

With the aim of quantifying the CO₂ output released to the atmosphere by TMCL, other authors have proposed an indirect method using continuous measurement of the partial pressure of CO₂ dissolved in the lake at a single point of observation with a NDIR CO₂ sensor⁴³. Figure 5 shows a comparison between the results obtained by the indirect method and the one showed in this work in the period 2013–2020. Both methods show a relatively good agreement with the exception of the period from August 2017 to the end of the time series: during 2017, the two methods showed an increase in the first half of the year that continued in the data reported in this work (further supported by the data measured at the flank of the volcano by the automatic geochemical station showed in Fig. 4), but is stopped drastically in the emission data obtained by the indirect method⁴³ (green line in Fg. 5). Two possible arguments or a combination of both can be made to explain this lack of correlation between both measurement methods in the second half of 2017: (1) the indirect method shows a systematic increase in the estimated CO₂ emission rate during the dry season of Taal (from November to June), and a systematic decrease in the rainy seasons, that might be caused by the so-called gas beracun^50,51; the process is driven by an stratification in acidic volcanic lakes caused by the entering of cold and fresh rainwater that coated the acidic lake waters with less saline, less acidic, and colder meteoric water. This process would lead a less efficient flushing of CO₂ through the lake’s top^50,51 and a worse or less accurate response of the indirect method to estimate the CO₂ emission rate from the entire surface of TMCL during the rainy season, from June to October; and (2) the loss of the correlation between the pCO₂ at the measuring site and the CO₂ emission rate of the TMCL after 2016, because, as was mentioned before, the measuring site might be located in an area that exhibits background emission values during the phase of unrest.

Figure 5

Temporal variations of the diffuse CO₂ degassing rate at the TMCL from 2013 to 2018 (solid red circles), the CO₂ flux estimated by the continuous measurement of the partial pressure of CO₂ dissolved in the lake at a single point of observation³⁹ (green line), and the trend of the seasonally adjusted CO₂ flux data reported by other authors³⁹ (blue line). Grey columns show the number of earthquakes per week³⁹.

The direct method to estimate the CO₂ emission rate from the TMCL used in this work, was able to detect a significant increase in 2017–2018 that represents the main and earliest precursor sign of the January 2020 eruption. The robustness of the 2017–2018 precursory increase of the diffuse CO₂ emission is supported also by the coeval increase observed in the automatic geochemical station (Fig. 4). The soil CO₂ efflux values measured at the automatic geochemical station suggests that the volatiles exsolved in the third magmatic intrusion likely started to reach the surface of the volcano on March 2017, but strain and stress changes in the crust due to such magma rise at depth were able to produce a significant increase of the seismic activity roughly 17 months after (Fig. 5). Other pulses of CO₂ emission cannot be ruled out in Taal due to subsequent magma movements due to the absence of diffuse CO₂ emission data after August 2019. The excellent agreement between the CO₂ emission values measured at TMCL and the soil CO₂ efflux values measured at the automatic geochemical station (Fig. 4), confirms that the latter methodology is an excellent complement to the gas emission surveys.

The gas emission data reported here, monitored with much higher detail since January 2016, have been very useful to provide a more reliable long-term warning of the 2020 eruption compared to the onset of precursory seismicity. The inspection of Fig. 3a suggests that between 2010 and 2012, magma rose to shallower crustal positions beneath Taal and that this upward magma movement was accompanied by a remarkable seismic unrest. Later, in 2014, degassing preceded the seismic unrest in few months. This observation might be due to a new slow and aseismic intrusion of magma in a shallower level that later exceeded the strength of the shallow host rocks⁵², or that the observed 2015 seismic unrest was triggered by the quiescent degassing caused by vesiculation or crystallization of the stagnant magma at shallower depths as suggested by other authors in different volcanic systems⁵³. During the third magmatic intrusion that led to the volcanic eruption, its initial phase started likely in March 2017 or before and was almost aseismic probably because it occurred at deeper levels and/or magma utilized existing open conduits produced by the two previous unrest stages. The process continued until November 2017, and later magma upward movements (until January 2020), were progressively causing the destabilization of the stagnant magma and opened the eruptive conduit, accompanied by the most energetic seismic activity of the period 2008–2020. Other geochemical parameters seem to support the occurrence of an injection of magmatic fluids into the Taal volcano hydrothermal system from the rise of fresh less and-degassed magma in 2017–2018⁵⁴. The temporal variations of the soil CO₂ efflux values measured at the automatic geochemical station preceded changes in the seismicity recorded by the Philippine Institute of Volcanology and Seismology (PHIVOLCS), similarly to what has been observed in other volcanic systems⁵⁵, and showed an excellent agreement with the discrete surveys performed in the same period (Fig. 4). This observation confirms the utility of the continuous monitoring of gas emission to monitor more accurately the timing of magma movements at depth.

Conclusions

The 10 years series of diffuse CO₂ degassing rate form TMCL, completed since 2016 with a continuous time series of soil CO₂ efflux in the north flank of the volcano, have been very useful to detect three early precursory signals of magmatic intrusion occurring beneath Taal volcano. Such magmatic rise led anomalous increase of the CO₂ emission rate from the TMCL measured in 2010–2011 and 2015 that preceded important changes in the seismic activity and other geophysical parameters of roughly 2–3 month, likely due to the decoupling between the amount of CO₂ exsolved from the ascending magma and the levels of the magmatic rise. New injection of magmatic volatiles into the hydrothermal system in 2017–2018 was caused by a third magmatic rise that likely occurred at deeper zones and could have been enough to trigger the 2020 volcanic eruption.

The geochemical data presented in this study represent the earliest warning precursor signals to the January 2020 eruptive activity. Both discrete CO₂ emission surveys and the continuous soil CO₂ efflux monitoring performed at a specific location have provided useful information to detect magma rise episodes. Continuous CO₂ flux monitoring definitely helps to accurately forecast volcanic eruptions, but regular CO₂ emission surveys should be promoted in order to complete the information as they cover wider areas of the volcano.

Methods

The statistical-graphical analysis of the data was based on representing the cumulative normal distribution of the data on a probability scale⁵⁶. On this scale, if a distribution is normal, the plot of the cumulative normal distribution versus the values results in a straight line. The direct reading of the variable at 50% probability provides the value of the average value ($\overline{{\text{x}}}$); the 84% probability reading provides the average value plus one standard deviation ($\overline{{\text{x}}} + \sigma$) and at 16% the average value minus one standard deviation ($\overline{{\text{x}}} - \sigma$). Similarly, $\overline{{\text{x}}} + 2\sigma$ and $\overline{{\text{x}}} - 2\sigma$ can be read at the 2% and 98% percentiles.

Diffuse CO₂ emission surveys from the TMCL (1.216 km²) were carried out performing an average 140 surface CO₂ efflux measurements per survey at the water–air interface and following the accumulation chamber method⁵⁷. The CO₂ efflux-meter was equipped with an non-dispersive IR CO₂sensor LICOR LI-820, able to measure in the range 0–2 mol%, with an accuracy of ~ 4%. The accumulation chamber connected to the non-dispersive IR CO₂ sensor was mounted on a floating device to allow the measurement at the water surface²⁵. The reproducibility of the CO₂ efflux-meter is 10% in the range 100−10,000 g m^–2 d^–1. This random error is based on the uncertainty calculated from the variability of the measurements carried out in the laboratory. In order to convert volumetric to mass flux rates, atmospheric pressure, temperature, and height of the accumulation chamber were taken into account. CO₂ efflux spatial distribution maps were constructed using conditional sequential Gaussian simulations. 100 equiprobable simulations were made by means of sGs algorithm^58,59, according to an experimental variogram model that fitted the experimental variogram. Each map is constructed by 3,041 square interpolated cells of 400 m² of surface and the CO₂ emission rate was estimated by the sum of the cells of the 100 simulations average map. The standard deviations of the 100 simulated values of total CO₂ output were assumed to be the characteristic values of its uncertainty⁵⁹. Continuous monitoring of soil CO₂ efflux measurements were performed by an automatic geochemical station installed in the northern sector of Taal volcano (14°1′14.3"N, 120°59′56.6"E) on 25 January 2016. A mechanical arm automatically placed a chamber over the ground every hour, and the soil CO₂ efflux was measured according to the accumulation chamber method⁵⁷ and an infrared sensor DRÄGER POLYTRON IR CO₂. Additionally, to avoid a possible influence of external parameters in the endogenous CO₂ emissions, soil water content and temperature at a depth of 40 cm and atmospheric parameters (wind speed and direction, air temperature and humidity, rainfall, and barometric pressure 1 m above the ground) were recorded simultaneously. All data were stored on an SD memory card and sent by GSM telemetry to the ITER-INVOLCAN facilities.