Seismic Data, Photographic Images and Physical Modeling of Volcanic Plumes as a Tool for Monitoring the Activity of Nevado del Ruiz Volcano, Colombia

Quantification of volcanic plume parameters is a fundamental task to characterize the behavior of an active volcano. The volcanic plume mass, flow rate and ash injection were determined from seismic data, in addition to photographic images and integration of scaling laws of several volcanic plume models, for the period from 1985 to 2017 for the Nevado del Ruiz Volcano (NRV), Colombia. With these parameters we quantified the ash volume emitted during this period and established a relationship between seismicity and the volcanic plume parameters. The results revealed a decrease of approximately two orders of magnitude in the volume of ash plumes from the November 13, 1985, eruption (0.12 km³) to the September 1, 1989, eruption (1.43 × 10⁻³ km³). This pattern continued for the June 30, 2012, eruption and 2015–2017 eruptive cycle, with volumes five times smaller than that observed in 1989. The results also exhibited a correlation between the radiated seismic energy (RSE) of the volcanic tremor and ash load for higher (>1 km) and longer-duration (>240 s) plumes. It was possible to calculate a minimum value of ash load based on RSE release and reduced displacement (RD, a means of normalizing volcanic tremors to a common scale) of volcanic tremor signals associated with the eruptions for the period 2015-2017. Moreover, changes in the volume of the ash plume were correlated with changes in the RD and RSE associated with different stages of volcanic activity. These findings can be used as a tool for monitoring the NRV. The continuously decreasing ash plume volumes from 1985 to 2017 suggest a common volcanic cycle that is almost ending. On the other hand, the evidence of new magmatic input in 2007 might suggest that a new volcanic cycle started on that date and is still in the process of ascending magma. It is likely that in the near future surface evidence of the new cycle will be observed at the NRV.

Introduction

Volcanic ash plumes pose a serious threat to both the neighboring population (e.g., Sparks et al., 1997; Horwell and Baxter, 2006; Jenkins et al., 2015) and civil aviation (e.g., Guffanti et al., 2010; Guffanti and Tupper, 2015). Their formation and dynamic behavior are controlled by the generation processes of ash (e.g., Dürig et al., 2012; Dellino et al., 2014; Dioguardi et al., 2016) and the fluid dynamic conditions at the source, which are described by the eruption source parameters (e.g., Wilson and Walker, 1987; Woods, 1988; Carazzo et al., 2008; Dellino et al., 2014). Quantification of volcanic plume parameters is fundamental to characterize the behavior of an active volcano and therefore a crucial step toward assessing its risk. Different approaches have been proposed. Starting from the seminal work of Sparks (1986), who proposed a basic physical model to calculate plume height, various 0D (0D models refer to scaling laws that relate plume height to eruptive mass or mass flow rate with no spatial variables being considered, e.g., Mastin et al., 2009), 1D (e.g., Devenish, 2013; Mastin, 2014), and 3D (e.g., Cerminara et al., 2016) models have been developed to quantify the main parameters of a volcanic plume (Bonadonna et al., 2002b; Kaminski et al., 2011; Degruyter and Bonadonna, 2012; Woodhouse et al., 2013; Suzuki and Koyaguchi, 2015; see Costa et al., 2016 for a review).

On the other hand, the current availability of cameras, videos, and other devices is becoming useful to understand in real-time the dynamics of volcanic plumes. In addition, seismic data can supply important information about the inner activity of a volcano. By combining all these data together, it is possible to obtain a better image of the volcano behavior, helping to monitor and forecast increases in volcanic activity.

Various studies have used such data to obtain a picture of the eruptive history of different volcanoes. For example, at Mt. Etna (Italy), Andronico et al. (2013) related ash emissions and the seismo-acoustic signals associated with them. They conclude that it is possible to differentiate two ash emissions types from the seismo-acoustic signals; one related with closed and the other with open conduits. McNutt et al. (2013) studied the 2009 explosive eruptions of Redoubt Volcano (Alaska) using seismic and infrasound signals. They associated the seismic energy and pressure obtained from infrasound data with column height. They found a correlation of seismic energy with gas release (SO₂) and pressure with column height. Dürig et al. (2015) analyzed the 2010 eruption of Eyjafjallajökull Volcano (Iceland) by using high-resolution video and aerial observations. They were able to observe the pulsatile activity of the volcano and to quantify the velocity of the pulses of ash emission. In addition, they estimated mass flux and plume height in good agreement with data from plume height and mass discharge models. De Angelis et al. (2016) studied the gas and ash explosions at Santiaguito Volcano (Guatemala) and the associated seismic signals (acoustic waveforms) and thermal infrared images with the aim to assess the bulk density of the eruptive plume, in addition to the fraction of ash and gas in the eruptive plume. They concluded that small to moderate explosions contained small fractions of ash. Romero et al. (2016) studied the eruption dynamics of Tungurahua volcano (Ecuador) based on fieldwork, thermal images and photographic images. They modeled the source parameters and suggested changes in the eruptive style. Fee et al. (2017) estimated the erupted mass of Sakurajima Volcano (Japan) by using infrasound waveforms, in combination with ash and gas data. They inverted infrasound waveforms associated with eruptions to quantify eruption flow rate and masses of 49 explosions. Their results agree with those obtained from ground-based ash collection and SO₂ data. Although far from complete, this brief list demonstrates the importance of using integrated approaches with multiple different data to assess key eruption source parameters.

Mastin et al. (2009) considered that one of the most exact methods to determine the duration of a volcanic eruption, even more precise than the direct observation, is the seismic signal associated with it. Such a seismic signal is called tremor, which is characterized by having a long duration and being sustained in time, with variations in seismic amplitude and frequency depending on the changes in the eruption dynamics (Denlinger and Moran, 2014). For this reason, in this study, the calculation of the duration of all the eruptions was based mainly on their associated seismic signals and checked against photographic records when available.

On the other hand, meteorological information such as atmospheric profiles has recently become a key factor to determine eruptive plume parameters more precisely (Woodhouse et al., 2013). A few decades ago, these factors were neglected or simplified. Therefore, to obtain more realistic source parameters, it is mandatory to use such information and the models that incorporate the information (Kaminski et al., 2011; Degruyter and Bonadonna, 2012). According to the method of Kaminski et al. (2011), the height of the eruptive plume for plinian (subplinian to ultraplinian) eruptions can be calculated by using a model that includes the following: the atmospheric conditions (atmospheric stratification), the magma temperature, the mass flow emitted, an entrainment coefficient that is a function of the plume buoyancy and environment temperature, and a partition factor related to the efficiency of magma fragmentation. For transient vulcanian eruptions, the Druitt et al. (2002) and Bonadonna et al. (2002a) models treat the plume as a discrete thermal and suggest that its height (H) can be obtained from the mass (M) emitted instantaneously into the atmosphere instead of the mass flow rate.

NRV is located in the center of Colombia and is well known for the deadly phreato-magmatic eruption on November 13, 1985. That eruption was cataloged as VEI = 3 based on some unclear estimations of duration and column height available at that time (Naranjo et al., 1986). On September 1, 1989, a VEI = 2 phreato-magmatic eruption occurred. On May 29 and June 30, 2012, two small phreato-magmatic eruptions were recorded. From 1985–1991 to 2015–2017, continuous small (vulcanian) phreatic eruptions occurred. Estimation of the mass flow rate, among other parameters, for most of those volcanic plumes has not been performed yet. In this study, we combine seismic records of the eruptions with ground-based photographic images and integral equations (0D models) or scaling laws derived from 1D plume models to estimate the source parameters for Nevado del Ruiz Volcano eruptions. With these results, we expect to increase the knowledge of the dynamics of the volcanic plumes of the NRV and contribute to volcano monitoring and volcanic hazard assessment. In addition, we expect a better quantification of plume parameters to enable comparisons between them and seismic parameters. Considerations regarding the 2010–2017 reactivation of the system are also presented.

Summary of the Recent Activity of Nevado Del Ruiz Volcano

The NRV is a stratovolcano located in center Colombia, with a summit height of approximately 5,311 m above sea level (a.s.l.). It is well known for the deadly and catastrophic phreatomagmatic eruption of November 13, 1985, which killed more than 20,000 people. In September 1989, another phreatomagmatic eruption occurred, with no victims. From 1986 to 1993, frequent small vulcanian eruptions were common at the NRV. In 2002 and 2003, an increase in seismicity and phreatic activity occurred. After almost 8 years of quiescence, it exhibited signs of reactivation in October 2010 (Londoño, 2016). In May and June 2012, two small phreatomagmatic eruptions occurred. After the end of 2014 and during 2015, continuous small vulcanian eruptions were common at the NRV. In 2015, a small dome was emplaced at the bottom of the active crater. This new activity is similar to that of the 1985–1991 period, with the difference of a dome building. Currently (Sept. 2018), the activity of the NRV continues at high levels.

New injection of magma in 2007 in the NRV zone was evidenced by changes in seismicity, deformation, and 3D seismic tomographic images (Londoño, 2016). The onset of the new activity was marked by increases in seismicity and SO₂ emission in October 2010, followed by a small ash emission. Furthermore, in 2011, a deep deformation source was detected far away to the SW of the volcano (Lundgren et al., 2015). In April 2012, strong shallow seismicity was detected to the SW of the active crater, and then on May 29 and June 30, two small phreatomagmatic eruptions occurred. The chemical composition of the ash was basically the same as the compositions of the ash from the 1985 to 1989 eruptions (Martinez et al., 2012). After those eruptions, the spatial distribution of volcano seismicity changed; several new volcano-tectonic (VT) seismogenic zones appeared, some of them located at the intersection of fault systems crossing the volcano. Between 2013 and 2014, some VT earthquakes were felt in Manizales City (30 km away from the NRV), some of them reaching local magnitudes (M_L) of between 4.7 and 5.0. In November 2014, small vulcanian eruptions characterized by ash emissions were common. Since July 2015, some changes in seismicity and deformation have been observed near the active crater in addition to a number of small vulcanian eruptions, which became more continuous. In September of the same year, a small dome was emplaced at the bottom of the active crater and is still growing (June 2018). After mid-2017, decreases in seismicity and the number of vulcanian eruptions compared to previous years (2015–2016) were observed. Recently, some authors have suggested that during from 2015 to 2016, new batches of fresh magma were injected into the plumbing system of the NRV (Londoño and Kumagai, 2018).

Methods

In this work, we used several scaling laws of volcanic plumes to estimate the source parameters for the activity of the NRV. For those eruptions with durations lasting < 4 min (240 s), we used the models of Druitt et al. (2002) and Bonadonna et al. (2002a), since they are established for an instantaneous release of mass (Bonadonna et al., 2002a). For eruptions longer than 4 min, we applied the models that are established for a steady plume of Kaminski et al. (2011), Degruyter and Bonadonna (2012) and Woodhouse et al. (2013).

In the Druitt et al. (2002) and Bonadonna et al. (2002a) models, the thermal plume height (H in m) is

$\begin{array}{cc}H=1.89J^4\text{ m }{(\phi Mc\left(T-To\right))}^{0.25},& (1)\end{array}$

where M is the mass of the solid material (kg), φ is the fraction of particles that contribute to the thermal mass of the plume (0.8), T-To is the difference in temperature between the atmosphere and the crater (To) and the emitted solid (T), and c is the specific heat of the solid (approximately 1,100 J kg⁻¹ K⁻¹). From Equation (1), it is possible to calculate H (in m) for a vulcanian plume as follows:

$\begin{array}{cc}\begin{array}{r}H=B{M}^{0.25}+Hv,\end{array}& (2)\end{array}$

where Hv is the height of the crater (in m) and B = 55 m kg⁻⁴ for Montserrat conditions (φ = 0.8; c = 1,204 J kg⁻¹ K⁻¹; T-To = 800 K). Since Montserrat volcano is similar in composition and eruptive style to the NRV, we used this formulation to model the NRV's vulcanian plumes. For the case of the NRV, an average of magma temperature of 1,173 K (Melson et al., 1990; Sigurdsson et al., 1990) and an average atmospheric temperature of 263 K based on meteorological data from IDEAM were used, given a value of B = 57 m J⁻¹.

In the model of Kaminski et al. (2011), which is valid for subplinian to ultraplinian eruptions, there is no wind considered, but there is a variable entrainment coefficient, whereas in the other models mentioned above, the entrainment coefficient is constant (0.1). In this model, the plume height (H in m) is

$\begin{array}{cc}\begin{array}{r}H=300\text{ m }{\text{s}}^4{\text{kg}}^{-4}*{Q_f}^{\frac{1}{4}},\end{array}& (3)\end{array}$

for H < 12,000 m and

$\begin{array}{cc}\begin{array}{r}H=5530\text{ m}+160\text{ m }{\text{s}}^4\mathrm{\text{k}}{g}^{-4}{Q_{\mathit{\text{f}}}}^{1/4},\end{array}& (4)\end{array}$

for H > 12,000 m,

where Q_f is the effective mass flow of the plume (ash + gas) (kg/s), which is given by

$\begin{array}{cc}\begin{array}{r}{Q}_f=\left[n_o+\phi (1-n_o)\right]Q_o\end{array}& (5)\end{array}$

where n_o is the gas fraction in magma, which ranges from 3 to 7% for silica magmas (Kaminski et al., 2011), and Q_o is the mass flow (kg/s), ϕ = f_p/100, where f_p is the partition factor, which depends on magma fragmentation (f_p = 100% for plinian eruptions and f_p <10% for effusive basaltic eruptions). For those eruptions in which portions of the emitted mass are not ejected into the atmosphere through the buoyant and ascending eruptive plume, but through pyroclastic flows or lavas, it is necessary to calculate the partition factor (see Equations 19–22 of Kaminski et al., 2011). From the effective mass flow (Q_f), it is possible to calculate the ash flow (Q_ash):

$\begin{array}{cc}\begin{array}{r}{Q}_{ash}=\phi (1-n_o)Q_f,\end{array}& (6)\end{array}$

where ϕ = f_p/100, and Q_f can be obtained from 1D models of eruptive plumes, such as

$\begin{array}{cc}\begin{array}{r}{Q}_f=a{H}^4+b,\end{array}& (7)\end{array}$

where a and b are coefficients that depend on atmospheric conditions. According to the works regarding plumes and turbulent jets by Carazzo et al. (2008) for tropical zones, such as that found in the NRV area, a = 7 × 10⁻¹¹ kg s⁻¹ m⁻⁴ and b = 0 kg s⁻¹ for H <18,000 m, and a = 2.78 × 10⁻¹¹ kg s⁻¹ m⁻⁴ and b = 2.5 × 10⁷ kg s⁻¹ for 18,000 <H <25,000 m. The partition factor, f_p, was taken to be 100% for all the eruptions of the NRV for the period from 1985 to 2017, except for the November 13, 1985, eruption, in which pyroclastic flows were generated (Calvache, 1990), corresponding to a f_p value of 98%. This is the only eruption with evidence of pyroclastic flows at the NRV in the studied period.

Kaminski et al. (2011) established a functional model for ash flows >1 × 10⁵ kg/s to determine the amount of ash at the top of the plume for a variable entrainment coefficient (see Figure 3 of Kaminski et al., 2011). Due to the possibility that for the NRV, several small phreatic eruptions had ash flow lower than 1 × 10⁵ kg/s, we obtained a fit to the curve in Figure 3 of Kaminski et al. (2011) with a power law, considering ash flows (Q_ash) lower than that value (E. Kaminski, pers. com. 2017). Consequently, the ash load (L given in mg m⁻³) for lower values of ash flows for the NRV was defined as

$\begin{array}{cc}\begin{array}{r}L=396.45\text{mg }{\text{m}}^{-3}\times ln(Q_{ash}/Q_o)-1,550.7{\text{mg m}}^{-3},\end{array}& (8)\end{array}$

where Q_o = 1 kg s⁻¹, a normalization constant.

An alternative model to calculate the mass flow of the volcanic plume is the model of Degruyter and Bonadonna (2012), which considers the local atmospheric conditions instead of a general atmospheric model, that is, it considers the wind conditions at the moment of the eruption. In that model, the mass flow is defined as

$\begin{array}{cc}\begin{array}{r}{Q}_o=\pi \frac{{\rho }_{a0}}{g^{\prime }}\left(\frac{2^{5/2}{\alpha }^2{\overline{N}}^3}{z_1^4}{H}^4+\frac{{\beta }^2{\overline{N}}^2\overline{v}}{6}{H}^3\right),\end{array}& (9)\end{array}$

where

$g^{\prime }=g\left(\frac{cT-c_{a0}{\theta }_{a0}}{c_{a0}{\theta }_{a0}}\right);\overline{N^2}=\frac{1}{H}{\displaystyle {\int }_o^H{N}^2\left(z\right)dz;}$

$\begin{array}{cc}\begin{array}{rcl}{N}^2(z)& =& \frac{g^2}{C_{a0}{\theta }_{a0}}\left(1+\frac{c_{a0}}{g}\frac{d{\theta }_a(z)}{dz}\right),\end{array}& (10)\end{array}$

$\begin{array}{cc}\overline{v}=\frac{1}{H}{\displaystyle {\int }_0^{H}v\left(z\right)dz},& (11)\end{array}$

where ρ_a0 is the reference density of the surrounding atmosphere, α is the radial entrainment coefficient (α = 0.1), β is the wind entrainment coefficient (β = 0.5), c_a0 is the heat capacity of the surrounding atmosphere, θ_a0 is the temperature of the surrounding atmosphere, g is the gravitational acceleration, g' depends on g and is computed via Equation (10), θ_a(z) is a profile of environment temperature, v(z) is a wind profile, $\overline{N}$ is the average buoyancy frequency of the atmosphere, $\overline{v}$ is the wind velocity across the plume height, z₁ is the maximum nondimensional height, and z is the vertical coordinate above the source.

A model similar to that of Degruyter and Bonadonna (2012) is the model of Woodhouse et al. (2013), in which the mass flow is given by

$\begin{array}{cc}\begin{array}{r}{Q}_o={\left[\frac{H}{0.318(\frac{1+1.373{\tilde{W}}_s}{1+4.266{\tilde{W}}_s+0.3527{{\tilde{W}}_s}^2})}\right]}^{3.95},\end{array}& (12)\end{array}$

where ${\tilde{W}}_s=1.44V_1/{\overline{N}H}_1$, V₁ is the wind velocity at the tropopause, H₁ is the local height of the tropopause (in m), and N is the atmospheric buoyancy frequency. Whereas, the model of Degruyter and Bonadonna (2012) considers the plume ascending in a calm atmosphere and bending over when it finds a strong wind field, the model of Woodhouse et al. (2013) considers the plume ascending in a linear shear crosswind in an intermediate regime, where the ascending plume speed and the wind speed are similar. For the NRV, we estimated the height of the tropopause and stratosphere as 12 and 20 km above the top of the volcano, respectively, from atmospheric models and radiosonde data from NOAA.

We calculated the density of the mixture contained in the volcanic plume (ρ_P) using the following formulation (Woods, 1988; Mastin, 2007):

$\begin{array}{cc}\begin{array}{r}{\rho }_P={\left[n_o{\rho }_g^{-1}+\frac{(1-n_0)}{{\rho }_m}\right]}^{-1},\end{array}& (13)\end{array}$

where ρ_m is the magma density (ash) and ρ_g is the gas density (vapor = 0.2 kg/m³ at 5,500 m a.s.l.).

In the previous 0D models, H refers to the height of the centerline of the plume, which is not identical to the height of the plume top in a bent-over situation. The 0D model of Mastin et al. (2009) established the column height based on an average eruption rate for different eruptions around the world. The column height (H in m) was defined as

$\begin{array}{cc}\begin{array}{r}H=2x10^3{V}^{0.241},\end{array}& (14)\end{array}$

where V is the volumetric flow rate (m³ of dense rock equivalent, DRE per second). In this model, both the column height at the top of the plume and column height at the center of the umbrella were considered (Mastin et al., 2009).

The radiated seismic energy (RSE) for the seismic signals (tremor) associated with the eruptions can be calculated by using the seismic power (Dibble, 1974; Cristofolini et al., 1987; Alparone et al., 2003):

$\begin{array}{cc}RSE=\pi {\rho }_s{v}_s{\left(r\text{ x }{10}^{-2}\right)}^2{A}^{2}t,& (15)\end{array}$

where ρ_s is the density of the upper part of the volcanic edifice (2,600 kg m⁻³, Londoño et al., 2014), v_s is the average seismic velocity of the shallowest layer (2,500 m s⁻¹, Londoño and Sudo, 2002), r is the source-station distance (in cm, assuming the active crater as the source), A is the average particle velocity (amplitude) filtered at the predominant frequency of the volcanic tremor, and t is the duration of the eruption.

Moreover, the instantaneous reduced displacement (RD, Fehler, 1983) for surface waves (Denlinger and Moran, 2014) for the volcanic tremor associated with the eruptions can be calculated using the following expression:

$\begin{array}{cc}\begin{array}{r}RD=\frac{A_{pp}\sqrt{\lambda r}}{2\sqrt{2}M},\end{array}& (16)\end{array}$

where A_pp is the maximum amplitude peak-to-peak (in cm) of the tremor filtered at the dominant frequency f (in Hz), λ is the wavelength for surface waves (in cm; λ = v_r/f, where v_r is the average velocity of surface waves = 2.16 km s⁻¹), and M is the magnification of the instrument. The factor 2√2 is the correction of the mean square root of the amplitude. Table 1 presents a list of the variables used in this study.

TABLE 1

Table 1. Description of the variables used in this study.

Data and Processing

The seismic, photographic, and SO₂ data used in this work belong to Servicio Geológico Colombiano (SGC). We used the seismic signals associated with the eruptions for the period from 1985 to 2017. Unfortunately, not all the eruptions before 2012 have photographic records, but most of them have data of direct visual observation obtained by Volcanological Observatory of Manizales from SGC. After 2012, most of them have photographic records. We used meteorological data regarding the studied zone belonging to IDEAM (the Colombia meteorological institute) and the NOAA Satellite and Information Service from USA. Additionally, we used data for some eruptive plumes height from Washington VAAC (Volcano Ash Advisory Code).

We used the maximum wind speed in the tropopause for the NRV by using the radiosonde located at Bogota City, 130 km to the SE of the crater. In a few cases, we used data from radiosondes located in Panamá, Ecuador or Curacao. We also calculated temperature gradients for the stratosphere of the NRV area.

For the period from 1985 to 1991, we analyzed only the larger eruptions, since there is partial information about column height and seismic records for many small phreatic eruptions that occurred. Nevertheless, we consider that the larger eruptions allow us to have an idea about the mass and emitted volumes between 1985 and 1991, which we estimated to account for approximately 85 to 90% of the total mass. We analyzed the volcanic plumes of the eruptions that occurred on September 11, 1985, November 13, 1985, January 4, 1986, May 4, 1986, July 20 and 29, 1986, June 10, 1987, March 25, 1988, September 1, 1989, and Apr 9, 26 and 29, 1991. Between 1992 and 2011, no eruptions were detected at the NRV. For the period from 2012 to 2017, a complete dataset is available, and almost all of the eruptions were analyzed.

Four hundred and twenty eruptive plumes were analyzed (see Appendix A). Eruptions that occurred during nighttime were not considered or included in this study; fortunately, these were few, and all of them were of small size and duration. Therefore, the results are not affected very much by this exclusion.

To measure the volcanic plume heights from the ground-based photographic images, we established a calibration with several places of known height observed in the images, taking the average height from available images of different photographic cameras for each eruption. Figure 1 shows an image of the location of the NRV and some examples of eruptive plumes. The error in the estimation of plume height was up to 10% (several 100 m); consequently, a higher error will occur for the parameters of the eruptive plume, such as the mass flow rate, mass, and volume. For each plume dataset, the source parameters were calculated using the different volcanic plume scaling laws mentioned above. Moreover, the values of some variables used in those models (Table 1) yield additional uncertainties of 10%; for instance, temperature can be affected by the presence of water for phreatic and phreato-magmatic eruptions. Thus, the total uncertainty of the source parameters calculated in this work can be estimated to be approximately 30% or even greater. Nevertheless, we are interested in the systematic evolution of the system as a function of time, which makes the absolute determination of the parameters less relevant.

FIGURE 1

Figure 1. (A) The location of the NRV and seismic (circles), camera (rectangles) and SO₂ (triangles) stations used for the analysis. (B) Examples of seismic records (digitized analog and digital) of some eruptions at the NRV (see Appendix C for details about the seismic signals). (C) Examples of typical eruptive plumes for the period from 1985 to 2017 (photographic images from the SGC database).

To measure the duration of the eruption, we used seismic records of the eruption. The duration of the eruption was constrained by both the seismic signal associated with the eruption and the synchronized photographic images available. The decay of the amplitude of the seismic signal was the main parameter used to define the end of the eruption, which in most of the cases was consistent with the stopping of ash emission observed in the photographic cameras. In a few cases, the end of the emission could not be observed directly from the cameras, but the seismic record was still available (see the example in Appendix B).

Additionally, the radiated seismic energy (RSE) and reduced displacement (RD) were calculated for the seismic signals associated with the eruptions. Figure 2 shows an example of a seismic signal associated with a small eruption and the parameters used to calculate the RD and RSE.

FIGURE 2

Figure 2. Example of parameters used for calculation of DR and RSE from seismic signals associated with eruptions. Amplitude values are in nm/s already corrected by instrumental response (name of variables are listed in Table 1).

Results

Plume Modeling

For the eruption of November 13, 1985, Naranjo et al. (1986) estimated a column height of 31 km, a duration of 20 min and a volume of dense material of 0.039 km³ based on modeling of the column height (Carey et al., 1986) and on eyewitness reports. We revisited that calculation and found that the data can be more precisely estimated. First, the duration of the eruption was corrected. Naranjo et al. (1986) estimated the duration of the eruption based on eyewitness observations, but the eruption occurred at nighttime, which makes the eyewitness technique not fully reliable. On the other hand, the signal recorded at several seismic stations around the volcano was used to calculate the duration of the eruption more accurately (Appendix C). By analyzing the recent high-resolution digitized analog seismic records available at SGC for the November 13, 1985, eruption carefully, we estimated that the eruption main pulse should have lasted between 90 and 105 min. After the main pulse, smaller pulses occurred and were recorded seismically. However, we only focus on the main phase of the eruption. Regarding the plume height, according to the work of Krueger et al. (1990) based on images of SO₂ from the Total Ozone Mapping Spectrometer (TOMS) instrument on the Nimbus 7 polar orbiting satellite available for the November 13 eruption, the estimated maximum column height reached up to 25 km and the main portion of the plume was between 7 and 14 km. If we remove the volcano height (approximately 5.3 km), we obtain a maximum eruptive plume height of approximately 20 km above the top of the volcano, which is at least 10 km less than the Naranjo et al. (1986) estimation, which was based on modeling. Krueger et al. (1990) used satellite images, which were possibly not available when Naranjo et al. (1986) and Carey et al. (1986) calculated the column height. Additionally, we estimated an atmospheric profile for the NRV region for the November 13, 1985, eruption based on data from NOAA from the Panamá radiosonde. We obtained a maximum wind speed at the troposphere for the Colombia, Panamá and Venezuela region (the region of plume dispersion) between 15 and 20 m/s and a temperature gradient of −2°C/km for the stratosphere. The heights of the tropopause and stratosphere were estimated to be 12 km and 20 km above the top of the volcano, respectively.

With these new data, the estimated flow rate ranged from 3.55 × 5 10⁷ ± 1.06 × 10⁷ kg/s according to the Mastin et al. (2009) scaling law to 6.69 × 10⁷ ± 2.01 × 10⁷ kg/s according to the Kaminski et al. (2011) scaling law. Therefore, the estimated mass emitted by the November 13, 1985, eruption was 3.6 × 10¹¹ ± 10¹¹ kg, 2.9 × 10¹¹ ± 8.9 × 10¹⁰ kg, 2.9 × 10¹¹ ± 8.7 × 10¹⁰ kg, and 5.7 × 10¹⁰ ± 1.7 × 10¹⁰ kg according to the models of Kaminski et al. (2011), Degruyter and Bonadonna (2012), Woodhouse et al. (2013), and Mastin et al. (2009), respectively.

For the September 1, 1985, eruption, we proceed in a similar manner. Méndez and Valencia (1991) calculated a volume of 1.6 × 10⁻³ km³ and a plume height of 6 to 8 km above the vent by using isopach data, in addition to a duration of 2 h 24 min. Based on analyzing in detail the digitized high-resolution analog seismic signal that was recently made available, the eruption lasted approximately 90 min for the main phase and approximately 6 h in total (Appendix C). We obtained an atmospheric profile for that day from radiosonde data from the BOGOTA station of the NOAA, and we computed a maximum wind speed in the troposphere of 9.7 km/s and a temperature gradient of −2°C/km for the stratosphere. The height of tropopause and stratosphere was estimated to be 12 and 20 km above the top of the volcano, respectively. We assumed a plume height of 8 km (Méndez and Valencia, 1991).

With the new data, the September 1, 1989, eruption emitted a mass of 2.7 × 10⁹ ± 8.2 × 10⁸ kg, 3.8 × 10⁹ ± 1.2 × 10⁹ kg, 5.2 × 10⁹ ± 1.5 × 10⁹ kg, and 4.3 × 10⁹ ± 1.3 × 10⁹ kg according to the models of Kaminski et al. (2011), Degruyter and Bonadonna (2012), Woodhouse et al. (2013), and Mastin et al. (2009), respectively.

For the minor eruptions, we proceeded in a similar manner, also. Figures 3, 4 show the mass flow rate (kg/s) and volume (m³), respectively, of eruptive plumes of the NRV for the period from 1985 to 2017 using different models. Figure 5 shows in detail those parameters for the period from 2015 to 2017. As it can be observed from Figure 3, the flow rate was variable during the studied period, with the highest values in 1985 and decreasing over time. A similar tendency was observed for the volume. For the period from 2015 to 2017 (Figure 5), those parameters exhibited a variable tendency, increasing and decreasing randomly. The total mass for the period from 2010 to 2017 for thermals (duration < 4 min) was 1.6 × 10⁸ ± 4.8 × 10⁷ kg, and for the other plumes, it ranged from 2.5 × 10⁸ kg to 5.39 × 10⁸ kg. The total plume volume was 1.8 × 10⁵ m³, using a mixture density of 3.3 kg/m³ (Ripepe et al., 2013), average air density of 0.7 kg/m³ at a height of 6 km, gas fraction between 3 and 6% for dacitic-andesitic rocks, and density of solid rock (without pores) between 2,400 and 2,600 kg/m³ for the NRV ash (Londoño et al., 2014; L. Martínez pers. com., 2017). Table 2 presents a summary of the source parameters obtained for the most relevant eruptions for the period from 1985 to 2017 for the NRV using different scaling laws.

FIGURE 3

Figure 3. The mass flow rate for eruptive plumes of the NRV with duration longer than 240 s and H>2 km for the period from 1985 to 2017. Each plume is represented by four overlapping bars. Y-axes in log scale. Error bars are included.

FIGURE 4

Figure 4. The volume (DRE) of the eruptive plumes of the NRV for the period from 1985 to 2017. The volume was calculated using the model of Degruyter and Bonadonna (2012) for eruptions with durations >240 s (red circles) and H>2 km (above the top of the vent) using the model of Druitt et al. (2002) for volcanic thermals with duration < 240 s (yellow circles). Y-axes in log scale. Error bars are included.

FIGURE 5

Figure 5. The mass flow rate (A) and volume DRE (B) of the eruptive plumes for the period from 2015 to 2017 for the NRV. The flow rate was calculated for eruptions with duration >240 s and H >2 km (above the vent). Y-axes in log scale. Error bars are included.

TABLE 2

Table 2. A comparison of eruptive source parameters for the NRV calculated with different volcanic plume models.

Additionally, we estimated the ash load (mg/m³) at the top of the eruptive plume (Kaminski et al., 2011) for the period from 1985 to 2017 using Equation (8). Figure 6 shows the results. In general, the ash load was bigger for the November 13, 1985, eruption (3,587 mg/m³). For the other eruptions, the ash load values ranged from approximately 500–3,400 mg/m³.

FIGURE 6

Figure 6. The ash load at the top of the eruptive plume for the eruptions of the NRV from 1985 to 2017. Error bars are included.

Relation Between Seismic Data and Eruption Plume Parameters

Figure 7 shows a plot of H vs. RD and H vs. RSE and a contour plot of H as a function of RD and RSE for the period from 2015 to 2017. In general, there was no correlation between H and RD or between H RSE for small values of H (< 3 km). In contrast, values of H >5 km corresponded to the highest values of RSE and RD, although they were not plotted. For small plumes (H < 3 km) with duration >240 s, there is a tendency of H as a function of RD and RSE (Figure 7C), that is, higher values of RD and RSE corresponded to high values of H, although there were several eruptions with H >1.5 km with relatively low RSE values.

FIGURE 7

Figure 7. The relationship between H (< 3 km) and RSE (A) and RD (B) for the studied eruptions from 2010 to 2017. Horizontal axes are in log scale. Vertical lines represent error bars. (C) Contour plot of H (< 3 km) and duration > 240 s as a function of RD and RSE.

Figures 8, 9 show details of the time series of volcanic plume parameters vs. RSE and RD, respectively, of tremors for the period from 2015 to 2017. From these figures, it is possible to observe that, in general, RSE and RD did not exhibit any clear relationship with plume parameters during that period. On the other hand, almost all the higher values of RD were registered during 2015, whereas higher values of RSE were distributed from 2015 to 2016.

FIGURE 8

Figure 8. Time series of radiated seismic energy (RSE) of volcanic tremor vs. eruptive plume parameters for the NRV for the period from 2015 to 2017. Different ranges of RSE are represented by colored squares. In the plume heights panel, all available data were plotted; for the rest of parameters, only eruptions with H >1 km and duration >240 s were plotted.

FIGURE 9

Figure 9. Time series of reduced displacement (RD) of volcanic tremor vs. volcanic plume parameters for the NRV for the period from 2015 to 2017. Different ranges of RD are represented by colored squares. In the plume heights panel, all available data were plotted; for the rest of parameters, only eruptions with H >1 km and duration>240 s were plotted.

We related the RSE with RD and the duration of the eruption with respect to the date of the eruption (Figure 10). There was a different pattern distribution of RSE and DR with respect to the duration of the eruption, depending on the date. On the other hand, the higher values of RD corresponded to the shortest durations of eruptions, which is not useful to forecast eruptions or volcanic behavior. This is a first indication that RSE is the preferable parameter when quantifying the volcanic tremor and its plume parameters compared to RD (see below).

FIGURE 10

Figure 10. The relationship between the duration of the eruption and the RSE (A) and DR (B) of the volcanic tremor associated with the eruption. The colors of squares represent different time periods (see Figure 12).

Discussion

According to Degruyter and Bonadonna (2012) and Woodhouse et al. (2013), if wind conditions are not accurately estimated or are neglected in the modeling of volcanic plumes, the mass flow rate can be underestimated by up to an order of magnitude. On the other hand, the duration of an eruption is another key parameter that must be determined as accurately as possible, since mass and volume parameters depend on it. Seismic records are a powerful tool to obtain eruption durations accurately (Mastin et al., 2009). In this study, we have used wind profiles and photographic and seismic records to calculate the eruptive plume parameters for the most recent period of activity of the NRV (1985–2017). It is noteworthy to mention that temperature is another factor that affects the eruption mass flow rate, as we pointed out previously. For phreatomagmatic pulses, it is possible that the magmatic temperature is less than that used in this study, due to presence of water; in addition, other factors, such as the initial thermal energy and mass of surface water, can be difficult to model for this type of eruption (Koyaguchi and Woods, 1996); therefore, we assume that there is another source of overestimation of such source parameters in our data not considered in this study.

With this in mind, the mass values obtained for the November 13, 1985, eruption imply that the volume of DRE was between 0.07 and 0.11 km³ when using a density of solid rock of 2,500 kg/m³ for the NRV (Melson et al., 1990). If we add the volume of the pyroclastic flows (0.009 km³, Calvache, 1990), the total volume of the November 13, 1985, eruption was at least 0.12 km³. This value is greater than that previously calculated by Naranjo et al. (1986), which was 0.03 km³ and that calculated by Calvache (1990), which was 0.02 km³. With the new data regarding the erupted volume, the November 13, 1985, eruption reaches a VEI of 4 (the lower limit of VEI = 4).

For the September 1, 1989, eruption, if we assume a density of solid rock of 2,500 kg/m³ (according to the compositional results of Méndez and Valencia, 1991), the volume of dense rock was 1.4 × 10⁻³ ± 1.0 × 10⁻⁴ according to the model of Degruyter and Bonadonna (2012). This volume agrees with that obtained by Méndez and Valencia (1991), who used isopach data.

Additionally, it is possible to establish an empirical relationship between RD and RSE with the ash load (L) at the top of the plume. If we neglect those volcanic plumes with H < 1 km, that is, the smallest size eruptions with low ash load values, we can construct a fit to obtain a minimum ash load (L_min) value (mg/m³), knowing the RD (in cm²) or RSE (in Joules) of the volcanic tremor signal associated with the eruptive column, as follows:

$\begin{array}{cc}\begin{array}{r}{L}_{min}=68.575{\text{mg}/\text{cm}}^3{\text{cm}}^{-2}\times \text{RD}{\text{cm}}^2+394.4{\text{mg}/\text{cm}}^3\end{array}& (17)\end{array}$

$\begin{array}{cc}\begin{array}{r}{L}_{min}=5\times 10^{-4}{\text{mg}/\text{cm}}^3{\text{cm}}^{-2}\times {\text{RSEcm}}^2+366.5{\text{mg}/\text{cm}}^3\end{array}& (18)\end{array}$

We choose to fit the minimum value instead of, for instance, an average, since we obtained a wide range of ash load values for similar RD or RSE. In this sense, the minimum value of ash load represents the lower limit of ash load we can obtain for a plume height, although it is possible to obtain higher values for the same plume height, up to the limit of theoretical values by Kaminski et al. (2011). Figure 11 shows this relationship and the fitted curves.

FIGURE 11

Figure 11. The relationship of ash load at the top of the eruptive plume with the RD (A) and RSE (B) of the associated tremor signal. The red circles represent the eruptive plumes with H >1 km. The red dotted lines represent the fitting curve for the minimum ash load value with respect to the RD and RSE.

As it can be observed from Figure 11, there is a minimum limit on the ash load depending on the value of RD or RSE; that is, the higher the RD or RSE value, the higher the minimum value of the ash load. Although it is possible to obtain a wide range of ash load values for the same value of RD or RSE, it appears that there is a lower limit on ash load for that value. For instance, for an RSE value of 1 × 10⁶ J, the minimum ash load value will always be greater than approximately 850 mg/m³. The same relation is valid for RD; for instance, a RD value of 6 cm² corresponds to a minimum value of approximately 750 mg/m³, and not less than that value. This finding is very interesting for volcano monitoring and forecasting the minimum ash load that an eruption of the NRV will contain based on the seismic signal only.

Moreover, there is a temporal variation of the RSE and RD of volcanic tremors related with the volume of ash for the period from 2015 to 2017. Figure 12 shows the comparison. From this figure, it is possible to observe three different stages or changes in RD and RSE with the cumulative volume of ash. The first stage (I in Figure 12) from March to the end of August 2015 exhibited a regular volume emission of ash, whereas RSE presented increases, and RD was steady, suggesting a semisealed magmatic system interacting with the hydrothermal system, partially blocking the output of solid material to the atmosphere. The second stage (II in Figure 12) from the end of August to the beginning of October 2015 exhibited an important increase in ash volume emission associated with a concurrent increase in RD and RSE, interpreted as a less sealed but pressurized magmatic system as a response of a dome emplacement at the crater bottom during September 2015 (SGC, 2015). The third stage (III in Figure 12) from March 2016 to the end of October 2016 exhibited an important increase in ash volume emission, whereas RD and RSE remained relatively steady, although RD exhibited a slight increase, implying a more open magmatic system, allowing solid material to be output freely, with a low amount of seismic energy needed to expel it. These stages seem to reasonably explain the current activity of the NRV (July 2018), which is characterized by low ash emissions associated with low energy seismic signal, whereas the dome is growing slowly, indicating that the volcanic system is almost open, allowing the ascent of solid material easily. These findings have some important implications for risk assessment: very strong eruptions, larger than the one of November 13, 1985, at the NRV, would probably require a drastic change in the conduits, such as blocking or pressurization. Such an imminent event should then be manifested in a significant increase in seismicity at the crater, in addition to deformation signals. On the other hand, the currently open condition of the volcanic system bears the possibility that volcanic eruptions of smaller or medium strength can occur without any considerable changes in seismicity or deformation due to the current open condition of the volcanic system. Therefore, it is mandatory to continue monitoring the activity of the NRV with uttermost vigilance and precision.

FIGURE 12

Figure 12. The time series showing the relation of the cumulative ash volume with the (A) reduced displacement (RD) and (B) radiated seismic energy (RSE) of the volcanic tremor associated with the eruption. Numbers I (March to end of August 2015), II (August to beginning of October 2015) and III (March to end of October 2016) correspond to three different stages of activity (see text for details).

On the other hand, the SO₂ flux is one of the most intriguing parameters observed at the NRV over time (Williams et al., 1990). With the aim to observe any relationship of source parameters and SO₂ release, we compared the mass flow and mass of eruptive plumes with those of SO₂ for the period 2015–2017, which was a period with continuous SO₂ measurements with DOAS instruments (Figure 1). Figure 13 shows this comparison. The SO₂ flux sometimes increased when the mass flow increased, but in other cases exhibited a contrary tendency. A decrease in SO₂ flux in December 2015 does not correspond to any eruption with H >2 km. In contrast, during May and June 2016, an increase in SO₂ was not associated with any eruption, likely suggesting passive degassing in the NRV during that time, that is, release of large amounts of SO₂ without an eruption.

FIGURE 13

Figure 13. A comparison of mass flow rate (red triangles) of volcanic plumes with duration >240 s and H >2 km calculated with the model of Mastin et al. (2009) with SO₂ flux (gray lines) for the NRV, 2015 to 2017. Each plotted SO₂ value corresponds to a measurement obtained every 5 min on average. Error bars are included.

Moreover, with the new calculated volume for the November 13, 1985, eruption, it is possible to revisit the question of the discrepancy between the SO₂ excess released to the atmosphere and the volume of the eruption. According to Krueger et al. (1990), the amount of SO₂ emitted by the eruption was 7x10⁸ kg; according to Williams et al. (1990), the magma volume needed to explain such SO₂ is approximately 0.92 km³, whereas Sigurdsson et al. (1990) estimated 0.3 km³. The new volume for this eruption, as mentioned previously, was 0.12 km³, approximately seven times less than expected according Williams et al. (1990) and only 2.5 times less than the value calculated by Sigurdsson et al. (1990). This discrepancy is less than the factor of 7–30 hitherto accepted and calculated by Naranjo et al. (1986) and Calvache (1990). Although the new volume still is less than that needed to explain all the SO₂ released, it is in the range of the magma bodies size beneath the NRV derived from seismic velocity anomalies from a 3D seismic tomography (Londoño and Sudo, 2002; Londoño and Kumagai, 2018). According to those seismic tomographic results, it is possible to estimate a shallow magma body (2–3 km depth) of approximately 10 km³ and another deep magma body (6–8 km depth) of approximately 50 km³ for 1985. In addition, to explain the constant degassing at the NRV, Williams et al. (1990) suggested a minimum magma volume between 4.6 and 9.2 km³ beneath the volcano, values that are in agreement with those calculated using seismic tomography for the shallow magma body. If we consider all the SO₂ released to the atmosphere by the NRV from 1985 to 2017, which was approximately 1.5 × 10¹⁰ kg, we can envision one or several magmatic reservoirs of approximately 60 km³ in total, in agreement with the estimated volume (summing both magmatic bodies) obtained from seismic tomography. With this new estimated volume, the ratio between the SO₂ emitted to the atmosphere and the ejected magmatic material is 5.1 × 10⁻⁴, which is two orders of magnitude less than that calculated by Williams et al. (1986). This ratio is similar to or smaller than those calculated for other volcanoes, such as St. Helens (USA) and El Chichón (Mexico) (Williams et al., 1986). On the other hand, if we compare the SO₂ released and the erupted magma for the November 13, 1985, eruption to other volcanoes according to the work of Wallace (2001), the new relation for the NRV fits better in the range of andesitic magma (Figure 14). Previously, these results corresponded to the basaltic magma range, which disagrees with the observed products in the field, classified as andesites and dacites (Melson et al., 1990; Sigurdsson et al., 1990).

FIGURE 14

Figure 14. The relation between SO₂ and the volume of erupted magma for several well-studied eruptions. (A) Agung, 1963; (EC) El Chichón, 1982; (Et) Etna, 1975; (F) Fuego, 1974; (L) La'scar, 1989; (Lk) Laki, 1783-1784; (Lq) Lonquimay, 1989; (ML) Mauna Loa, 1984; (M) Mount St. Helens, 1980; (Pc) Pacaya, 1972; (P) Pinatubo, 1991; (Rb) Rabaul, 1994; (R) Redoubt, 1989-90; (Rz) Ruiz, 1985; (S) Stromboli; and (Bt) Bishop Tuff. The red circle (Rz_new) corresponds to the results of this work. Modified and adapted from Wallace (2001).

It has been argued that the November 13, 1985, eruption of the NRV was modest (Giggenbach et al., 1990; Krueger et al., 1990; Williams et al., 1990), but it is possible that it was not as small as it has been assumed. Several facts support this idea. First, the meteorological conditions that day were not favorable to remotely observe the eruptive plume by weather satellites, and there are not available satellite images for the NRV region at the moment of the eruption according to a search for GOES-6, METEOSAT-2 or GMS-3 satellite images using the National Centers for Environmental Information NOAA web browser.

(https://www.ncdc.noaa.gov/has/HAS.FileAppRouter?datasetname=3645&subqueryby=STATION&applname=&outdest=FILE); Second, SO₂ plumes were observed remotely 15 h later (Krueger et al., 1990) by dedicated satellites that can detect particles of fine ash not observed by ground-based radar and eyes. The SO₂ plume was detected more than 1,000 km away from the volcano to the Atlantic Ocean, suggesting a coexisting very large amount of probably fine ash ejected to the atmosphere, but not deposited close to the volcano. Third, three decades after the eruption we obtained information about a volcanic ash layer of that eruption at least one or two mm thick as far as 215 km far from the volcano (geologist Italo Reyes eyewitness and photographic record, SGC photographic database) in the Belencito municipality in Boyacá Department (120 km to NE of Bogotá), suggesting that the isopachs were probably much more extended than previously calculated by Naranjo et al. (1986). Fourth, the volume of pyroclastic current density (PCD) products could not be accurately estimated, since much of that material was incorporated into the lahar, as pointed out by Calvache (1990). In addition to these facts, the reevaluation of the duration of the eruption performed in this work leads us to conclude that the eruption was at least one order of magnitude larger than previously determined. It is possible that much of the ejected material was incorporated into the lahars (coarse material) or deposited far away (fine material) up to the Atlantic Ocean, it being impossible to realistically estimate the true ejected volume of magmatic material. The newly available data help to constrain the real size of this eruption with a higher precision.

The estimated volume in this work for the September 1, 1989, eruption is in agreement with that obtained by Méndez and Valencia (1991). For this eruption, more data were available, including geological (petrographic, stratigraphic), geophysical (seismic) and observational data, supporting the idea that the combination of different datasets leads to more realistic and consistent volcanic plume parameters.

The estimated volume for the small eruptions of May 29 and June 30, 2012 exhibited a discrepancy with that calculated by Martinez et al. (2012) using an isopachs approach. For the May 29 and June 30, 2012 eruptions, the ash volumes calculated by Martinez et al. (2012) were 1.59 × 10⁶ m³ and 5.83 × 10⁴ m³, respectively, whereas in this study, we calculated a volume between 7.6 × 10⁴ (with the model of Mastin et al., 2009) and 1.4 × 10⁵ m³ (with the model of Woodhouse et al., 2013) for May 29 and between 2.2 × 10⁵ m³ (with the model of Mastin et al., 2009) and 3.1 × 10⁵ m³ (with the model of Woodhouse et al., 2013) for the June 30, 2012 eruptions. This difference may be due to different facts. First, the estimated plume height for these eruptions could be affected by a bent-over situation due to wind at the moment of the eruption; therefore, it is possible that the volume estimated can be biased by the uncertainty in the plume height. For the May 29 eruption, the wind speed at 10 km a.s.l. was 18 m/s, whereas for June 30, it was 24 m/s, according to IDEAM meteorological reports. With this source of error for column height in mind, according to the model of Mastin et al. (2009), the volume of the May 29 eruption was half an order of magnitude less than that of 30 June 2012. It is possible that the volume calculated using the isopachs approach, which uses the plume height as the umbrella-cloud height, was underestimated for the June 30 eruption. Based on the available data and reports, the volcanic plume height of the June 30 eruption was higher than that of May 29; according to the VAAC homepage (Volcanic Ash Advisory, 2010; http://www.ssd.noaa.gov/VAAC/ARCH12/RUIZ/2012E291325.html) and a report from an airplane pilot, the plume height of the May 29 eruption (3:10 am local time) was 5.7 km above the top of the volcano. According to photographic images of SGC, the column height of the 30 June eruption (5:47 pm local time) was at least 7 km above the top of the volcano. Second, the seismic signal of the June 30 eruption was longer than that of May 29. Based on these facts, we conclude that the June 30 eruption was larger than the May 29 eruption. The differences in tephra dispersion may be explained by two facts: first, heavy rain was falling all night long on June 30 in a wide region of the Caldas, Risaralda and Quindío Departments (central Colombia), probably removing part of the fine ash fall in several places, leading to underestimate the tephra dispersion for that eruption (Martinez et al., 2012). Second, it is possible that during the May 29 eruption, the conduit was more sealed than June 30, which involved more pressure and gas, ejecting fine ash over larger distances than the June 30 eruption. The difference in the seismic amplitudes of the signals supports this conjecture; although the June 30 eruption long lasted more than May 29, the latter exhibited a higher instantaneous seismic amplitude (higher RD) than the former (Appendix B).

A continuous decrease in the volume of ash plumes of NRV of about two orders of magnitude from the November 13, 1985 eruption to the September 1, 1989 eruption, and of about 5 orders of magnitude from 1989 to 2012, suggests that a volcanic cycle is ending. Moreover, the correlations between RSE and RD of the volcanic tremor and the ash load, as well as a correlation between RSE, RD and the ash volume for the studied period, allow us to divide the activity of NRV in different stages. These findings could be used to better understand the behavior of NRV in the future.

Our study shows that a holistic analysis of datasets from multiple different sources, such as records of gas, seismicity, and local meteorological and atmospheric conditions, in combination with photographic images and scaling laws based on physical plume models, leads to a significantly improved assessment of eruption source parameters, hence being a promising tool for volcanic monitoring and risk assessment, not only for NRV but also to other active volcanoes.

Concluding Remarks

Volcanic plume parameters such as the mass flow rate, erupted mass, ash flow rate, and volume were calculated for the most recent eruptive period of the NRV, 1985–2017, by integrating different datasets and using different integral volcanic plume models. These new data suggest that the eruption was larger than previously determined. With these new data, the November 13, 1985, eruption had a VEI = 4, with a volume of 0.12 km³. The eruption lasted longer than initially assumed, based on a detailed analysis of seismic signals associated with the eruption. Moreover, it is possible to establish a relation between the RSE of volcanic tremors and the ash load at the top of the eruptive plume, which allows the RSE of volcanic tremors to be used to provide advice about the risk of volcanic ash for aviation, helping in the volcanic risk assessment of this region.

In addition, the ongoing reactivation of the NRV started in October 2010 and corresponds to the same eruptive cycle of 1985. It is likely that the same magmatic body is acting and releasing large amounts of SO₂ to the atmosphere, in addition to volcanic ash, which has been gradually decreasing over time, suggesting that this cycle is ending. On the other hand, it is argued that currently, a new deep magmatic body is intruding to the S of the NRV (15–18 km, Lundgren et al., 2015), which probably started to move up in 2007 (Londoño, 2016) and was recently highlighted at the surface by the emplacement of a small dome at the bottom of the active crater. This phenomenon indicates that the NRV still is a very active volcano, with the possibility of new eruptions in the near future. Knowing of the current behavior of a volcano by means of the dynamics of the volcanic plume and its relationship with the seismicity is a useful approach to elucidate future behavior of that volcano. In this study that approach was used for NRV and hopefully it can be applied to other volcanoes.

Author Contributions

JL had the main idea, wrote portions of the manuscript and processed SO₂ and plume parameters. BG processed the seismic data, collected height plumes from digital images and wrote some portions of the manuscript.

Funding

This work was partially supported by Universidad Católica de Manizales, Servicio Geológico Colombiano, and Universidad de Caldas, Manizales.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We want to thanks to C. Bonadonna and W. Degruyter for their help and for providing computer codes for calculation of volcanic plume parameters. We also thank C. Bonadonna for her careful review of and comments on this paper, E. Kaminski for his suggestions and comments about his model and application to small volcanic plumes, and our colleagues at Servicio Geológico Colombiano (SGC) at the Volcanological and Seismological Observatory of Manizales for their support and suggestions. This work was partially supported by the SGC Direction of Geohazard, Universidad de Caldas and Universidad Católica de Manizales. Luca de Siena, E. Kaminski, and another reviewer, in addition to the associate editor, Yosuke Aoki, and the chief editor, Valerio Acocella, made important suggestions and comments that considerably improved the final manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2018.00162/full#supplementary-material

References

Alparone, S., Andronico, D., Lodato, L., and Sgroi, T. (2003). Relationship between tremor and volcanic activity during the Southeast Crater eruption on Mount Etna in early 2000. J. Geophys. Res. 108:2241. doi: 10.1029/2002JB001866

CrossRef Full Text | Google Scholar

Andronico, D., Castro, L., Sciotto, M., and Spina, L. (2013). The 2010 ash emissions at the summitcraters of Mt Etna: relationship with seismo-acous-tic signals. J. Geophys. Res. 118, 51–70. doi: 10.1029/2012JB009895

CrossRef Full Text | Google Scholar

Bonadonna, C., Macedonio, G., and Sparks, R. S. J. (2002a). “Numerical modelling of tephra fallout associated with dome collapses and Vulcanian explosions: application to hazard assessment on Montserrat,” in The Eruption of Soufriere Hills Volcano, Montserrat, from 1995 to 1999, eds T. H.Druitt and B. P. Kokelaar (London: Memoirs, Geological Society) 21, 517.2–537.2.

Google Scholar

Bonadonna, C., Mayberry, G. C., Calder, E. S., Sparks, R. S. J., Choux, C., Jackson, P., et al. (2002b). “Tephra fallout in the eruption of Soufrière Hills Volcano, Montserrat,” in The Eruption of Soufriere Hills Volcano, Montserrat, from 1995 to 1999, eds T. H.Druitt and B. P. Kokelaar (London: Memoirs) 21, 483–516.

Google Scholar

Calvache, M. L. (1990). Geology and Volcanology of the Recent Evolution of Galeras Volcano, Colombia. M.S. Thesis, Louisiana State University.

Carazzo, G., Kaminski, E., and Tait, S. (2008). On the dynamics of volcanic columns: a comparison of field data with a new model of negatively buoyant jets. J. Volcanol. Geotherm. Res. 178, 94–103. doi: 10.1016/j.jvolgeores.2008.01.002

CrossRef Full Text | Google Scholar

Carey, S., Sigurdsson, H., Fritz, W., and Naranjo, J. L. (1986). The November 13, 1985 eruption of Nevado del Ruiz: tephra fall deposit and column height modelling (abstr.). EOS Trans. Am. Geophys. Union 67, 403.

Google Scholar

Cerminara, M., Esposti Ongaro, T., and Berselli, L. C. (2016). ASHEE-1.0: a compressible, equilibrium–Eulerian model for volcanic ash plumes. Geosci. Model Dev. 9, 697–730. doi: 10.5194/gmd-9-697-2016

CrossRef Full Text | Google Scholar

Costa, A., Suzuki, Y. J., Cerminara, M., Devenish, B. J., Ongaro, T. E., Herzog, M., et al. (2016). Results of the eruptive column model inter-comparison study. J. Volcano. Geotherm. Res. 326, 2–25. doi: 10.1016/j.jvolgeores.2016.01.017

CrossRef Full Text | Google Scholar

Cristofolini, R., Gresta, S., Imposa, S., Menza, S., and Patani, G. (1987). An approach to problems on energy sources at Mt. Etna.from seismological and volcanological data. Bull. Volcanol. 70, 729–736.

Google Scholar

De Angelis, S., Lamb, O. D., Lamur, A., Hornby, A. J., Aulock, F. W., Chigna, G., et al. (2016). Characterization of moderate ash-and-gas explosions at Santiaguito volcano, Guatemala, from infrasound waveform inversion and thermal infrared measurements. Geophys. Res. Lett. 43, 6220–6227. doi: 10.1002/2016GL069098

PubMed Abstract | CrossRef Full Text | Google Scholar

Degruyter, W., and Bonadonna, C. (2012). Improving on mass flow rate estimates of volcanic eruptions Geophys. Res. Lett. 39:L16308. doi: 10.1029/2012GL052566

CrossRef Full Text | Google Scholar

Dellino, P., Dioguardi, F., Mele, D., D'Addabbo, M., Zimanowski, B., Büttner, R., et al. (2014). Volcanic jets, plumes and collapsing fountains: evidence from large-scale experiments, with particular emphasis on the entrainment rate. Bull. Volcanol. 76:834. doi: 10.1007/s00445-014-0834-6

CrossRef Full Text | Google Scholar

Denlinger, R. P., and Moran, S. C. (2014). Volcanic tremor masks its seismogenic source: results from a study of noneruptive tremor recorded at Mount St. Helens, Washington. J. Geophys. Res. Solid Earth 119, 2230–2251. doi: 10.1002/2013JB010698

CrossRef Full Text | Google Scholar

Devenish, B. J. (2013). Using simple plume models to refine the source mass flux of volcanic eruptions according to atmospheric conditions. J. Volcanol. Geotherm. Res. 256, 118–127. doi: 10.1016/j.jvolgeores.2013.02.015

CrossRef Full Text | Google Scholar

Dibble, R. R. (1974). “Volcanic seismology and accompanying activity of Ruapehu Volcano, New Zealand,” in Physical Volcanology, eds L. Civetta, P. Gasparini, G. Luongo, and A. Rapolla (Amsterdam: Elsevier) 49–85. doi: 10.1016/B978-0-444-41141-9.50009-7

CrossRef Full Text | Google Scholar

Dioguardi, F., Dürig, T., Engwell, S. L., Gudmundsson, M. T., and Loughlin, S. C. (2016). “Chapter 6, Investigating source conditions and controlling parameters of explosive eruptions: some experimental-observational- modelling case studies,” in Updates in Volcanology - From Volcano Modelling to Volcano Geology, ed K. Nemeth (Zagreb: InTech), 141–162. doi: 10.5772/63422

CrossRef Full Text | Google Scholar

Druitt, T. H., Young, S. R., Baptie, B., Bonadonna, C., Calder, E. S., Clarke, A. B., et al. (2002). Episodes of cyclic Vulcanian explosive activity with fountain collapse at Soufrière Hills Volcano, Montserrat. Geol. Soc. London Memoir 21, 281–306. doi: 10.1144/GSL.MEM.2002.021.01.13

CrossRef Full Text | Google Scholar

Dürig, T., Gudmundsson, M. T., Karmann, S., Zimanowski, B., Dellino, P., Rietze, M., et al. (2015). Mass eruption rates in pulsating eruptions estimated from video analysis of the gas thrust-buoyancy transition—a case study of the 2010 eruption of Eyjafjallajökull, Iceland. Earth Planets Space 67:180. doi: 10.1186/s40623-015-0351-7

CrossRef Full Text | Google Scholar

Dürig, T., Sonder, I., Zimanowski, B., Beyrichen, H., and Büttner, R. (2012). Generation of volcanic ash by basaltic volcanism, J. Geophys. Res. 117:B01204. doi: 10.1029/2011JB008628

CrossRef Full Text | Google Scholar

Fee, D., Izbekovb, P., Kim, K., Yokoo, A., Lopez, T., Pratae, F., et al. (2017). Eruption mass estimation using infrasound waveform inversion and ash and gas measurements: evaluation at Sakurajima Volcano, Japan. Earth Planet. Sci. Lett. 480, 42–52. doi: 10.1016/j.epsl.2017.09.043

CrossRef Full Text | Google Scholar

Fehler, M. C. (1983). Observation of volcanic tremor at Mount, St. Helens volcano. J. Geophys. Res. 88, 3476–3484.

Google Scholar

Giggenbach, W. F. N., Garcia, P. A., Londofio, C. L., Rodriguez, V. N., Rojas, G., and Calvache, M. L. (1990). The chemistry of fumarolic vapor and thermal-spring discharges from Nevado del Ruiz volcanic-magmatic-hydrothermal system, Columbia. J. Volcan. Geotherm. Res. 42, 13–39. doi: 10.1016/0377-0273(90)90067-P

CrossRef Full Text | Google Scholar

Gil, F. (1991). Determinación de la Respuesta Instrumental de la Red Sismológica en el Volcán Nevado del Ruiz. INGEOMINAS. Internal Report. 94.

Guffanti, M., Casadevall, T. J., and Budding, K. E. (2010). Encounters of Aircraft With Volcanic Ash Clouds: A Compilation of Known Incidents, 1953-2009. Reston, VA: U.S. Department of the Interior; U.S. Geological Survey.

Guffanti, M., and Tupper, A. (2015). “Volcanic Ash Hazards and Aviation Risk,” in Volcanic Hazards, Risks and Disasters, eds J. F. Shroder and P. Papale (Boston, MA: Elsevier) 87–108. doi: 10.1016/B978-0-12-396453-3.00004-6

CrossRef Full Text | Google Scholar

Horwell, C. J., and Baxter, P. J. (2006). The respiratory health hazards of volcanic ash: a review for volcanic risk mitigation. Bull. Volcanol. 69, 1–24. doi: 10.1007/s00445-006-0052-y

CrossRef Full Text | Google Scholar

Jenkins, S. F., Wilson, T., Magill, C., Miller, V., Stewart, C., Blong, R., et al. (2015). “Volcanic ash fall hazard and risk” in Global Hazards and Risk, eds S. C. Loughlin S. parks, S. K. Brown S. F. Jenkins, and C. Vye-Brown (Cambridge, UK: Cambridge University Press), ISBN:978-1-107-11175-2. doi: 10.1017/CBO9781316276273.005

CrossRef Full Text | Google Scholar

Kaminski, E., Tait, S., Ferrucci, F., Martet, M., Hirn, B., and Husson, P. (2011). Estimation of ash injection in the atmosphere by basaltic volcanic plumes: The case of the Eyjafjallajokull 2010 eruption. J. Geophys. Res. 116:B00C02. doi: 10.1029/2011JB008297

CrossRef Full Text | Google Scholar

Koyaguchi, T., and Woods, A. W. (1996). On the formation of eruption columns following explosive mixing magma and surface-water. J. Geophys. Res. 101, 5561–5574. doi: 10.1029/95JB01687

CrossRef Full Text | Google Scholar

Krueger, A. J., Walter, L. S., Schnetzler, C. C., and Doiron, S. D. (1990). TOM's measurement of the sulfur dioxide emitted during the 1985 Nevado del Ruiz eruptions. J. Volcanol. Geotherm. Res. 41, 7–15.

Google Scholar

Londoño, J. M. (2016). Evidence of recent deep magmatic activity at Cerro Bravo-Cerro Machín volcanic complex, central Colombia. Implications for future volcanic activity at Nevado del Ruiz, Cerro Machín and other volcanoes. J. Volcanol. Geotherm. Res. 324, 156–168. doi: 10.1016/j.jvolgeores.2016.06.003

CrossRef Full Text | Google Scholar

Londoño, J. M., and Kumagai, H. (2018). 4D seismic tomography for Nevado del Ruiz volcano, Colombia, 2000–2016. J. Volcanol. Geotherm. Res. 358, 105–123. doi: 10.1016/j.jvolgeores.2018.02.015

CrossRef Full Text | Google Scholar

Londoño, J. M., Ponce, P., and Ortega, A. (2014). Un modelo Geofísico del Volcán Nevado del Ruiz. In Spanish. Internal Report. Servicio Geológico Colombiano, 1–35.

Londoño, J. M., and Sudo, Y. (2002). Velocity structure and a seismic model for Nevado del Ruiz Volcano (Colombia). J. Volcanol. Geotherm. Res. 119, 61–87. doi: 10.1016/S0377-0273(02)00306-2

CrossRef Full Text | Google Scholar

Lundgren, P., Samsonov, S. V., López, C. M., and Ordoñez, M. (2015). Deep source model for Nevado del Ruiz Volcano, Colombia, constrained by interferometric synthetic aperture radar observations. Geophys. Res. Lett. 42, 4816–4823. doi: 10.1002/2015GL063858

CrossRef Full Text | Google Scholar

Martinez, L. M., Zuluaga, I., and Ceballos, J. A. (2012). Informe Sobre la Actividad Eruptiva del Volcán Nevado del Ruiz: Mayo 29 y Junio 30 de 2012. Servicio Geológico Colombiano. OVS-Manizales. Informe interno. 67.

Mastin, L. G. (2007). A user-friendly one-dimensional model for wet volcanic plumes, Geochem. Geophys. Geosyst. 8, 1-24. doi: 10.1029/2006GC001455

CrossRef Full Text | Google Scholar

Mastin, L. G. (2014). Testing the accuracy of a 1-D volcanic plumemodel in estimating mass eruption rate. J. Geophys. Res. Atmos. 119, 2474–2495. doi: 10.1002/2013JD020604

CrossRef Full Text | Google Scholar

Mastin, L. G., Guffanti, M., Servranckx, R., Webley, P., Barsotti, S., Dean, K., et al. (2009). A multidisciplinary effort to assign realistic source parameters to models of volcanic ash-cloud transport and dispersion during eruptions. J. Volcanol. Geothermal Res. 186, 10–21. doi: 10.1016/j.jvolgeores.2009.01.008

CrossRef Full Text | Google Scholar

McNutt, S. R., Thompson, G., West, M. E., Fee, D., Stihler, S., and Clark, E. (2013). Local seismic and infrasound observations of the 2009 explosive eruptions of Redoubt Volcano, Alaska. J. Volcanol. Geotherm. Res. 259, 63–76. doi: 10.1016/j.jvolgeores.2013.03.016

CrossRef Full Text | Google Scholar

Melson, W. G., Allan, J. F., Jerez, D. R., Nelen, J., Calvache, M. L., Williams, S. N., et al. (1990). Water contents, temperatures and diversity of the magmas of the catastrophic eruption of Nevado del Ruiz, Colombia, November 13, (1985). J. Volcanol. Geotherm. Res. 41, 97–126. doi: 10.1016/0377-0273(90)90085-T

CrossRef Full Text | Google Scholar

Méndez, R. A., and Valencia, H. T. (1991). La Erupción del volcán Nevado del Ruiz el primero de septiembre de 1989: datos geológicos y modelo de la erupción. Boletín Geológico. Ingeominas, Bogotá 35, 22–42.

Google Scholar

Naranjo, J. L., Sigurdsson, H., Carey, S. N., and Fritz, W. (1986). Eruption of the Nevado del Ruiz Volcano, Colombia, on 13 November 1985: tephra fall and lahars. Science 233, 961–963. doi: 10.1126/science.233.4767.961

PubMed Abstract | CrossRef Full Text | Google Scholar

Ripepe, M., Bonadonna, C., Folch, A., Delle Donne, D., Lacanna, G., and Voight, B. (2013). Ash-plume dynamics and eruption source parameters by infrasound and thermal imagery: the 2010 Eyjafjallajokull eruption. Earth Planet Sci. Lett. 366, 112–121. doi: 10.1016/j.epsl.2013.02.005

CrossRef Full Text | Google Scholar

Romero, J. E., Douillet, G. A., Vargas, S. V., Bustillos, J., Troncoso, L., Díaz, J., et al. (2016). Dynamics and style transition of a moderate, Vulcanian-driven eruption at Tungurahua (Ecuador) on February 2014: pyroclastic deposits and hazard considerations. Solid Earth 8, 697–719. doi: 10.5194/se-2016-159

CrossRef Full Text | Google Scholar

SGC (2015). Informe de actividad volcánica Segmento Norte, Octubre 2015. In Spanish, Internal Report. Servicio Geológico Colombiano. 1–93. Available online at: https://www.sgc.gov.co/Publicaciones/Informes%20tcnicos/Informe%20del%20mes%20de%20octubre%202015.pdf

Sigurdsson, H., Cornell, W., and Carey, S. (1990). Influence of magma withdrawal on compositional gradients during the AD 79 Vesuvius eruption. Nature 345, 519–521. doi: 10.1038/345519a0

CrossRef Full Text | Google Scholar

Sparks, R. S. J. (1986). The dimensions and dynamics of volcanic eruption columns. Bull. Volcanol. 48, 3–15. doi: 10.1007/BF01073509

CrossRef Full Text | Google Scholar

Sparks, R. S. J., Bursik, M. I., Carey, S. N., Gilbert, J. S., Glaze, L., Sigurdsson, H., et al. (1997). Volcanic Plumes. John Wiley & Sons, Inc.

Google Scholar

Suzuki, Y. J., and Koyaguchi, T. (2015). Effects of wind on entrainment efficiency in volcanic plumes, J. Geophys. Res. Solid Earth 120, 6122–6140. doi: 10.1002/2015JB012208

CrossRef Full Text | Google Scholar

Volcanic Ash Advisory (2010) Volcanic Ash Advisory from London – Issued Graphics. VAAC (Washington Volcanic Ash Advisory Centre). Available online at: http://www.ssd.noaa.gov/VAAC/ARCH12/RUIZ/2012E291325.html.,

Wallace, P. J. (2001). Volcanic SO₂ emissions and the abundance and distribution of exsolved gas in magma bodies. J. Volc. Geotherm. Res. 108, 85–106. doi: 10.1016/S0377-0273(00)00279-1

CrossRef Full Text | Google Scholar

Williams, S. N., Stoiber, R. E., Garcia, N. P., Londono, A. C., Gemmell, B. J., Lowe, D. R., et al. (1986). Eruption of the Nevado del Ruiz volcano, Colombia, on 13 November, 1985: gas flux and fluid geochemistry. Science 233, 964–967. doi: 10.1126/science.233.4767.964

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, S. N., Sturchio, N. C., Calvache, M. L., Londono, A., and Garcia, N. (1990). Sulfur dioxide from Nevado del Ruiz volcano, Colombia: Geotherm. Res. 42, 53–68. doi: 10.1016/0377-0273(90)90069-R

CrossRef Full Text | Google Scholar

Wilson, L., and Walker, G. P. L. (1987). Explosive volcanic eruptions, VI. Ejecta dispersal in plinian eruptions: the control of eruption conditions and atmospheric properties. Geophys. J. Int. 89, 657–679.

Google Scholar

Woodhouse, M. J., Hogg, A. J., Phillips, J. C., and Sparks, R. S. J. (2013). Interaction between volcanic plumes and wind during the 2010 Eyjafjallajökull eruption, Iceland. J. Geophys. Res. 116, 92–109. doi: 10.1029/2012JB009592

CrossRef Full Text | Google Scholar

Woods, A. W. (1988). The fluid dynamics and thermodynamics of eruption columns. Bull. Volcanol. 50, 169–193. doi: 10.1007/BF01079681

CrossRef Full Text | Google Scholar