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 know the behavior of an active volcano. Volcanic plume mass, flow rate and ash injection were determined from integration of scaling laws of several volcanic plume models, seismic data, and photographic images for the period 1985–2017 for Nevado del Ruiz Volcano (NRV), Colombia, with the aim to quantify the ash emitted volume during this period and to establish a relationship between seismicity and those volcanic plume parameters. The results show a decrease of ~ two orders of magnitude in the volume of ash plumes from the November 13, 1985 eruption (0.12 km3) to the September 1, 1989 eruption (1.43x10-3 km3). This pattern is preserved for the June 30, 2012 eruption and 2015-2017 eruptive cycle, with volumes five times smaller than in 1989. The results also show a correlation between 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 radiated seismic energy (RSE) release and reduced displacement (RD) of volcanic tremor signal associated with the eruption, which could be used as a monitoring tool. Moreover, for the period 2015-2017 changes in volume of ash plume were correlated with changes in reduced displacement (RD, a way of normalizing volcanic tremor to a common scale) and RSE, associated with different stages of volcanic activity. The continuous decreasing of ash plume volumes from 1985 to 2017 suggests a same volcanic cycle 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 still in process of magma ascending. Probably in the next future surface evidence of the new cycle can be observed at NRV or other volcanoes around it.

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 3 ) to the September 1, 1989, eruption (1.43 × 10 −3 km 3 ). 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 realtime 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 2 ) 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 groundbased ash collection and SO 2 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 . 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 1985From -1991From to 2015From -2017 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 groundbased 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 2 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 andBonadonna (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 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 −1 K −1 ). From Equation (1), it is possible to calculate H (in m) for a vulcanian plume as follows: where Hv is the height of the crater (in m) and B = 55 m kg −4 for Montserrat conditions (ϕ = 0.8; c = 1,204 J kg −1 K −1 ; 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 −1 .
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 for H < 12,000 m and for H > 12,000 m, where Q f is the effective mass flow of the plume (ash + gas) (kg/s), which is given by 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 ): where φ =f p /100, and Q f can be obtained from 1D models of eruptive plumes, such as 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 −11 kg s −1 m −4 and b=0 kg s −1 for H<18,000 m, and a=2.78 × 10 −11 kg s −1 m −4 and b = 2.5 × 10 7 kg s −1 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 5 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 5 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 −3 ) for lower values of ash flows for the NRV was defined as 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 where 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, N is the average buoyancy frequency of the atmosphere, v is the wind velocity across the plume height, z 1 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 whereW s = 1.44V 1 /NH 1 , V 1 is the wind velocity at the tropopause, H 1 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): where ρ m is the magma density (ash) and ρ g is the gas density (vapor = 0.2 kg/m 3 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 where V is the volumetric flow rate (m 3 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): where ρ s is the density of the upper part of the volcanic edifice (2,600 kg m −3 , Londoño et al., 2014), v s is the average seismic velocity of the shallowest layer (2,500 m s −1 , 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: 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 −1 ), 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.

DATA AND PROCESSING
The seismic, photographic, and SO 2 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, 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 1992and 2011 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 groundbased 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.
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.

Plume Modeling
For the eruption of  estimated a column height of 31 km, a duration of 20 min and a volume of dense material of 0.039 km 3 based on modeling  ) 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   Table 1). 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 2 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 7 ± 1.06 × 10 7 kg/s according to the Mastin et al. (2009) scaling law to 6.69 × 10 7 ± 2.01 × 10 7 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 11 ± 10 11 kg, 2.9 × 10 11 ± 8.9 × 10 10 kg, 2.9 × 10 11 ± 8.7 × 10 10 kg, and 5.7 × 10 10 ± 1.7 × 10 10 kg according to the models of Kaminski   For the September 1, 1985, eruption, we proceed in a similar manner. Méndez and Valencia (1991) calculated a volume of 1.6 × 10 −3 km 3 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 9 ± 8.2 × 10 8 kg, 3.8 × 10 9 ± 1.2 × 10 9 kg, 5.2 × 10 9 ± 1.5 × 10 9 kg, and 4.3 × 10 9 ± 1.3 × 10 9 kg according to the models of Kaminski  For the minor eruptions, we proceeded in a similar manner, also. Figures 3, 4 show the mass flow rate (kg/s) and volume (m 3 ), 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 8 ± 4.8 × 10 7 kg, and for the other plumes, it ranged from 2.5 × 10 8 kg to 5.39 × 10 8 kg. The total plume volume was 1.8 × 10 5 m 3 , using a mixture density of 3.3 kg/m 3 (Ripepe et al., 2013), average air density of 0.7 kg/m 3 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 3 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.
Additionally, we estimated the ash load (mg/m 3 ) 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 3 ). For the other eruptions, the ash load values ranged from approximately 500-3,400 mg/m 3 .  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.

Relation Between Seismic Data and Eruption Plume Parameters
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).

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 . 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.
For the September 1, 1989, eruption, if we assume a density of solid rock of 2,500 kg/m 3 (according to the compositional results of Méndez and Valencia, 1991), the volume of dense rock was 1.4 × 10 −3 ± 1.0 × 10 −4 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 3 ), knowing the RD (in cm 2 ) or RSE (in Joules) of the volcanic tremor signal associated with the eruptive column, as follows: L min = 68.575mg/cm 3 cm −2 ×RDcm 2 + 394.4mg/cm 3 (17) L min = 5×10 −4 mg/cm 3 cm −2 ×RSEcm 2 + 366.5mg/cm 3 (18) 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.
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 6 J, the minimum ash load value will always be greater than approximately 850 mg/m 3 . The same relation is valid for RD; for instance, a RD value of 6 cm 2 corresponds to a minimum value of approximately 750 mg/m 3 , 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.
On the other hand, the SO 2 flux is one of the most intriguing parameters observed at the NRV over time . With the aim to observe any relationship of source parameters and SO 2 release, we compared the mass flow and mass of eruptive plumes with those of SO 2 for the period 2015-2017, which was a period with continuous SO 2 measurements with DOAS instruments (Figure 1). Figure 13 shows this comparison. The SO 2 flux sometimes increased when the mass flow increased, but in other cases exhibited a contrary tendency. A decrease in SO 2 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 2 was not associated with any eruption, likely suggesting passive degassing in the NRV during that time, that is, release of large amounts of SO 2 without an eruption.
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 2 excess released to the atmosphere and the volume of the eruption. According to Krueger et al. (1990), the amount of SO 2 emitted by the eruption was 7x10 8 kg; according to Williams et al. (1990), the magma volume needed to explain such SO 2 is approximately 0.92 km 3 , whereas Sigurdsson et al. (1990) estimated 0.3 km 3 . The new volume for this eruption, as mentioned previously, was 0.12 km 3 , 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 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).
all the SO 2 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 3 and another deep magma body (6-8 km depth) of approximately 50 km 3 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 3 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 2 released to the atmosphere by the NRV from 1985 to 2017, which was approximately 1.5 × 10 10 kg, we can envision one or several magmatic reservoirs of approximately 60 km 3 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 2 emitted to the atmosphere and the ejected magmatic material is 5.1 × 10 −4 , 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 2 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).
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 2 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 2 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 6 m 3 and 5.83 × 10 4 m 3 , respectively, whereas in this study, we calculated a volume between 7.6 × 10 4 (with the model of Mastin et al., 2009) and 1.4 × 10 5 m 3 (with the model of Woodhouse et al., 2013) for May 29 and between 2.2 × 10 5 m 3 (with the model of Mastin et al., 2009) and 3.1 × 10 5 m 3 (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 3 . 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 2 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 2 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.

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.