Mass Eruption Rates of Tephra Plumes During the 2011–2015 Lava Fountain Paroxysms at Mt. Etna From Doppler Radar Retrievals

Real-time estimation of eruptive source parameters during explosive volcanic eruptions is a major challenge in terms of hazard evaluation and risk assessment as these inputs are essential for tephra dispersal models to forecast the impact of ash plumes and tephra deposits. Between 2011 and 2015, Etna volcano has produced 49 paroxysms characterized by lava fountains generating tephra plumes that reached up to 15 km a.s.l.. We analyzed these paroxysms using the 23.5 cm wavelength Doppler radar (VOLDORAD 2B) signals along with visible camera images of the monitoring network of the Istituto Nazionale di Geofisica e Vulcanologia, Osservatorio Etneo. Range gating of the radar beam allows the identification of the active summit craters in real-time, no matter the meteorological conditions. The radar echoes help to mark (i) the onset of the paroxysm when unstable lava fountains, progressively taking over Strombolian activity, continuously supply the developing tephra plume, then (ii) the transition to stable fountains (climax), and (iii) the end of the climax with a waning phase, therefore providing paroxysm durations. We developed a new methodology to retrieve in real-time a Mass Eruption Rate (MER) proxy from the radar echo power and maximum Doppler velocity measured near the emission source. The increase in MER proxies is found to precede by several minutes the time variations of plume heights inferred from visible and X-Band radar imagery. A calibration of the MER proxy against ascent models based on observed plume heights leads to radar-derived climax MER from 2.96 × 104 to 3.26 × 106 kg s-1. The total erupted mass (TEM) of tephra was computed by integrating over beam volumes and paroxysm duration, allowing quantitative comparisons of the relative amounts of emitted tephra among the different paroxysms. When the climactic phase can be identified, it is found to frequently release 76% of the TEM. Calibrated TEMs are found to be larger than those retrieved by satellite and X-band radar observations, deposit analyses, ground-based infrared imagery or dispersion modeling. The radar-derived mass load parameters therefore represent a very powerful all-weather tool for the quantitative monitoring and real-time hazard assessment of tephra plumes at Etna.

Real-time estimation of eruptive source parameters during explosive volcanic eruptions is a major challenge in terms of hazard evaluation and risk assessment as these inputs are essential for tephra dispersal models to forecast the impact of ash plumes and tephra deposits. In this aim, taking advantage of the 23.5 cm wavelength Doppler radar (VOLDORAD 2B) monitoring Etna volcano, we analyzed 47 paroxysms produced between 2011 and 2015, characterized by lava fountains generating tephra plumes that reached up to 15 km a.s.l. Range gating of the radar beam allows the identification of the active summit craters in real-time, no matter the meteorological conditions. The radar echoes help to mark (i) the onset of the paroxysm when unstable lava fountains, taking over Strombolian activity, continuously supply the developing tephra plume, then (ii) the transition to stable fountains (climax), and (iii) the end of the climax, therefore providing paroxysm durations. We developed a new methodology to retrieve in real-time a Mass Eruption Rate (MER) proxy from the radar echo power and maximum Doppler velocity measured near the emission source. The increase in MER proxies is found to precede by several minutes the time variations of plume heights inferred from visible and X-Band radar imagery. A calibration of the MER proxy against ascent models based on observed plume heights leads to radar-derived climax MER from 2.96 × 10 4 to 3.26 × 10 6 kg s −1 . The Total Erupted Mass (TEM) of tephra was computed by integrating over beam volumes and paroxysm duration, allowing quantitative comparisons of the relative amounts of emitted tephra among the different paroxysms. When the climactic phase can be identified, it is found to frequently release 76% of the TEM. Calibrated TEMs are found to be larger than those retrieved by satellite and X-band radar observations, deposit analyses, ground-based infrared imagery, or dispersion modeling. Our methodology, potentially applicable to every Doppler radar, provides mass load parameters that represent a powerful all-weather tool for the quantitative monitoring and real-time hazard assessment of tephra plumes at Etna or any other volcano with radar monitoring.
Keywords: Etna, paroxysmal activity, Lava fountains, Doppler radar, mass eruption rate, total erupted mass INTRODUCTION Quantifying the so-called eruptive source parameters (Bonadonna et al., 2015) of tephra plumes is critical for hazard assessment of explosive volcanic eruptions and associated risk mitigation, as well as for a better understanding of the dynamics of eruption columns and plumes. The different eruptive source parameters are: the location of the eruptive vent, the start time and duration of an eruption, the plume height, the Mass Eruption Rate (MER), and the Total Grain Size Distribution (TGSD). These parameters are used by the Volcanic Ash Advisory Centers (VAACs) to initialize Volcanic Ash Dispersion and Transportation Models (VATDMs) in near real-time. A particularly challenging objective is to measure the MER in real-time. It is generally derived from empirical relationships between observed top heights of strong plumes and corresponding MERs inferred from scaling laws (Wilson et al., 1978;Sparks et al., 1997;Mastin et al., 2009). However, Mastin et al. (2009) and Degruyter and Bonadonna (2012) have reported that such empirical relationships between plume heights and MERs are subject to high uncertainties. MERs estimated from post-eruption deposits analyses themselves hold uncertainties highly dependent upon the selected methodology (Andronico et al., 2014a;Bonadonna et al., 2015;Spanu et al., 2016).
A way to operationally retrieve, i.e., in (near) real-time, the eruptive source parameters is to use remote sensing techniques. Radars represent particularly robust tools for realtime assessment of source parameters owing to their relatively high spatial resolution and acquisition rate, their all-weather detection capacity near the emission source allowing early warning and quantification. Fixed-beam transportable Doppler radars with high time resolution were for instance used to monitor and study the dynamics of Strombolian and mild Vulcanian activity, using either 23.5-cm wavelength radars mostly sensitive to lapilli-and block-sized tephras (Dubosclard et al., 1999(Dubosclard et al., , 2004Donnadieu, 2010, 2011;Donnadieu et al., 2011;Valade et al., 2012), or 1-cm wavelength micro rain radars well suited for lapilli and coarse ash detection (Seyfried and Hort, 1999;Hort et al., 2003;Scharff et al., 2015;Hort and Scharff, 2016). Strong Vulcanian to Plinian eruptions have also been surveyed with 5-cm (Harris and Rose, 1983) and 3-cm wavelength scanning weather radars (Marzano et al., 2013;Maki et al., 2016;Vulpiani et al., 2016). Those radars have shown their capabilities and strength to study the dynamics of tephra plumes in real time and to provide estimates of (some of) the source parameters a posteriori, although generally with a lack of output parameters crossvalidation.
At Etna (Figure 1A), one of the most active European volcanoes, the repetitive explosive activity and the risks associated with tephra plumes has led the Istituto Nazionale di Geofisica e Vulcanologia, Osservatorio Etneo (INGV-OE) to improve its monitoring network to better anticipate and measure Etna's ash emissions . The network is based on the use of different remote sensing measurements and VATDMs runs daily using fixed eruptive scenarios (e.g., Scollo et al., 2009Scollo et al., , 2010. In this context, a 23.5 cm-wavelength Doppler radar (VOLDORAD 2B) has been integrated into the INGV-OE monitoring network since 2009 (Donnadieu, 2012;Donnadieu et al., 2016). This radar recorded 43 out of 45 paroxysmal episodes from the New Southeast Crater (NSEC) between 2011 and 2013, and 4 from the Voragine Crater (VOR) in December 2015, totaling 47 paroxysms between January 2011 and December 2015 ( Figure 1B). Two paroxysms were missed on 19 July 2011 and 20 February 2013 due to power outage and radar maintenance. Paroxysms at Etna are powerful events lasting several hours and characterized by lava fountains generating high eruption columns accompanied or not by the emission of lava flows (Andronico et al., 2014a;Corsaro et al., 2017). The plumes typically reach 9-15 km above sea level, produce downwind fallout of lapilli (and sometimes bombs) up to several kilometers from the vents and ash fallout up to 400 km away from the volcano (Andronico et al., 2015). Considering that the mild Strombolian activity preceding the paroxysmal activity may last a few hours to a few days, and that the transition to a sustained tephra plume is currently not accurately predictable, it is crucial to quantify the evolution of the source parameters in near real-time.
Yet, measuring the whole set of eruptive source parameters of an eruption is not trivial. Indeed, in addition to the aforementioned observables, the Total Grain Size Distribution is also required (Bonadonna et al., 2015). The latter parameter is often incompletely estimated due to the limitations of tephra sampling near the summit craters. Indeed, the aforementioned paroxysms have two distinct fallout contributions. On the one hand, lava fountains are composed of dense ballistics and windpushed lighter blocks and lapilli that fall close to the source (i.e., <5 km). Despite the fact that they likely represent the dominant part of the total erupted mass, they are rarely sampled because the deposits are hardly distinguishable from those of previous eruptions and owing to recurrent fallout in the hardly accessible Valle del Bove (Andronico et al., 2014a;Spanu et al., 2016). On the other hand, lapilli and ash constituting the developing tephra plumes are often wind-drifted toward Southeast above the Ionian Sea, again preventing sampling. Incomplete deposit sampling leads, in turn, to high uncertainties on the retrieved Total Erupted Masses (TEMs), from which, the mean MERs are derived (Andronico et al., 2014a).
In this paper, we first describe the VOLDORAD 2B monitoring system and utilize the Doppler radar retrievals to qualitatively describe common features of the eruption dynamics during paroxysmal episodes of Etna between 2011 and 2015. We then present a new methodology to compute a proxy for the erupted mass only from the measured radar parameters. This methodology has potentially powerful application in realtime monitoring at Etna, but also at any volcano monitored by Doppler radars. Calibration of the mass proxy with plume ascent models parameterized with observed plume heights and with results from other methods leads to MER (potentially in realtime) and TEM estimates. Results are then discussed in the last section.
FIGURE 1 | (A) Mount Etna location and photograph of the summit crater (courtesy of Boris Behncke). Geometry of radar beam above Etna's summit craters: probed volumes are drawn at −3 dB, i.e., at half the power in the beam axis, and dashed lines indicates beam limit at −10 dB; (B) top view (after Oct. 10 2012); (C) S-N cross-section view (aperture angle of 8.3 • in elevation at −3 dB): before December 16 2013, 11 range gates (3,135-4,635 m) were monitored and 13 gates (2,685-4,485 m) after this date. Inset: for range gates above the emission source, the positive (v + max ) and negative (v − max ) radial velocities measured along-beam mainly stem from ascending and falling tephra, respectively.

MATERIALS AND METHODS
The VOLDORAD 2B Monitoring System Dubosclard et al. (1999Dubosclard et al. ( , 2004 carried out pioneering surveys at Etna using a transportable L-band fixed-beam radar showing the high potential of radars to quantitatively monitor explosive activity near the emission source. They found in particular a correlation between tremor amplitude and echo power and ejecta initial velocities. In 2009, a similar 23.5-cm wavelength radar, VOLDORAD 2B, was set up at La Montagnola station on Etna (2,610 m a.s.l.) with its fixed beam pointing to a zone right above the summit craters 3 km northward (Figures 1A,C). Since then, it has been continuously monitoring the tephra emissions in volumes close to the summit craters (Donnadieu, 2012;Donnadieu et al., 2015Donnadieu et al., , 2016. The 23.5 cm wavelength is well suited for the detection of lapilli and blocks/bombs allowing to probe inside the tephra column regardless of weather conditions. The high sampling rate (about 5 Hz) allows the realtime quantification and provides insight into the dynamics of the eruption column at time scales of individual explosions to that of entire eruptive/inter-eruptive periods. The radar beam is divided into successively probed 150 m-deep volumes (range gates) extending 1.2 km above the summit craters area along the N-S direction of the beam. This range gating provides spatial information on the explosive activity, allowing for instance the identification of the active crater or craters during simultaneous activity (Donnadieu, 2012). From 2009 to October 10, 2012, the radar beam aimed above the summit vent with azimuth and elevation (θ ) angles of 347.5 • and 13 • , respectively. After this date, the radar antenna was rotated to about 355.2 • in azimuth and 14.9 • in elevation ( Figure 1B) to better record the activity of the NSEC. On December 16 2013, two more proximal range gates were added, passing the number of recorded volumes ranging from 11 (3,135-4,635 m) to 13 (2,835-4,635 m; Figures 1B,C). VOLDORAD 2B simultaneously records the amplitude of the echo power backscattered by the tephra and their radial velocity (measured along-beam using the Doppler effect) in each range gate. Displays of radar parameters, the power spectral distribution as a function of radial velocities, are called Doppler spectra. Velocity component toward the radar are negative, and positive away from it (Sauvageot, 1992). For the range gates located above the emission source, the power associated with positive and negative radial velocities mainly stem from ascending and falling tephra respectively ( Figure 1C). Out of the time series of power and velocity parameters retrieved from the Doppler spectra (e.g., Dubosclard et al., 2004), two are most useful to quantify the mass loading of explosive activity, as explored in a following section: (i) the total power P(t) backscattered by tephra in each probed volume, which is directly related to the quantity and size of particles crossing the radar beam; and (ii) the maximum positive Doppler velocity v + max (t) as it can be geometrically related to the ejection velocities V(t) assuming vertical jets (Dubosclard et al., 1999(Dubosclard et al., , 2004Donnadieu, 2012;Scharff et al., 2015): with θ the elevation angle of the radar beam ( Figure 1). As θ was changed on 10 October 2012 from about 13 • to 14.9 • , V(t) ≈ 4.45v + max (t) for the 2011-2012 paroxysms, and V(t) ≈ 3.89v + max (t) afterwards, including the 2013-2015 paroxysms.

Plume Top Height Measurements
In order to retrieve absolute MERs from the radar parameters, we have used independently the MER obtained from the column height observations. In fact, the link between plume heights and mass eruption rates is one of the most studied among volcanic source parameters relationships (Mastin et al., 2009). Scollo et al. (2014) proposed a methodology to retrieve column heights at Etna from image analyses of the ECV visible camera (in Catania, 27 km away from Etna's summit craters), with an error of ±500 m . The method limitations include night and bad weather conditions preventing the use of this visible camera, and the maximum altitude of 9 km above sea level. In this case, the ECV measurements may be supplemented by satellite imageries to retrieve the maximum column height using the Dark Pixel procedure that assumes a thermal equilibrium between the plume top and the atmosphere (Wen and Rose, 1994;Prata and Grant, 2001;Corradini et al., 2016). When available, we have also used DPX4 X-band weatherradar data of the Italian civil protection, in addition to other remote sensor data (i.e., satellite and visible imagery) estimating the plume heights during the 23 November 2013 NSEC paroxysm and the December 2015 VOR Crater paroxysms (Corradini et al., 2016;Vulpiani et al., 2016).

The Radar Mass Eruption Rate Proxy
Several recent works using scanning weather radars aimed at estimating mass loading parameters of explosive eruptions. Marzano et al. (2006) produced a procedure to retrieve ash mass load parameters (i.e., VARR model) using an electromagnetic scattering model and Dual-polarization radar observables. Their work was applied to Etna paroxysmal activity in 2013 (Corradini et al., 2016;Montopoli, 2016) and in December 2015 (Vulpiani et al., 2016) using the volume information of the X-Band (3 cm wavelength) weather radar located at Catania airport (30 km south from the Etna's summit), with a 3-D scan time resolution of 10 min.
Taking advantage of the higher time (<0.1 s) and spatial (120 m) resolution of a fixed-beam radar similar to VOLDORAD 2B pointing right above the emission source, Gouhier and Donnadieu (2008) developed an inversion method based on the Mie Scattering Theory to retrieve the ejecta mass of individual outbursts during Strombolian activity at Etna in 2001. Because of their short emission time, Strombolian explosions were treated as quasi-instantaneous releases of particles in which all ejecta could be captured in the large volumes of the fixed beam during the recorded peak of echo power.
The continuous monitoring of Etna with the VOLDORAD 2B radar at high space-time resolution (150 m, 0.2 s) offers a good opportunity to estimate the mass load parameters of Etna paroxysms. However, the lack of accurate physical characterization of proximal tephra (i.e., from the lava fountaining) in terms of shape, size and density weakens assumptions on inputs to scattering simulations, in particular the particle size distribution, and brings out large uncertainty in the mass load outputs. Therefore, in the following, we present a new approach based on a simple analytical model to compute the tephra mass loading parameters from a mass proxy directly retrieved from Doppler radar observables above the vent, and then calibrated against values measured by other methods. Interestingly, this methodology does not require an accurate particle size distribution and is applicable to the most frequent cases of eruptions in which the tephra emission duration is longer than the time needed for tephra to cross the beam. It also has obvious application to improve real-time monitoring and hazard assessment of tephra plumes.
As our goal is to calibrate a (relative) mass proxy directly related to radar observables, the physical model does not need to mimic the complexity of the particle dynamics during the eruption but only to correlate with the MER evolution. In our simplified eruption model, spherical particles with a unique diameter D (in m), constant with time, cross the beam vertically at velocity V(t) (in m s −1 ) assumed equal to the maximum ejection velocity, and constant over the beam crossing height (Figure 2).
The number of ascending particles dN inside the volume probed during the radar sampling period dt between two successive measurements is therefore defined by: where n(t) is the number of particles per unit volume (m −3 ), and S is the entering surface area (m 2 ) of the jet into the beam, no matter its shape.
FIGURE 2 | Sketch of the physical model to compute the mass eruption rate from the power and velocity parameters measured by the radar. S defines the (arbitrary) entry surface of the ejecta at the bottom of the radar beam. V(t) is the ejected tephra velocity in the beam, assumed vertical, and v + max its component along-beam, as measured by the radar. θ is the radar beam elevation angle. P t (t) is the peak power transmitted into the atmosphere by the radar, and P + (t) the power backscattered from ejecta having a positive radial velocity, supposedly ascending vertically, in a range gate above the emission source.
Assuming spheres of density ρ (kg m −3 ), the total particle mass over time M (kg) is: Under the Rayleigh assumptions (D < λ/4 in Gouhier and Donnadieu, 2008; where λ is the 23.5 cm-wavelength of our radar), the power (P + (t), in mW) backscattered by ascending spheres homogeneously distributed inside a probed beam volume above the emission source can be obtained from the radar equation (e.g., Sauvageot, 1992): γ being a constant (mW m −3 ) gathering known parameters specific to the radar. The combination of Equations (1-4) leads to the timeintegrated mass of tephra M (kg) expelled through the probed volume between times t 1 and t 2 : where M * is the above integral and C (kg mW −1 m −1 ) the constant factor before it. Under the complete Mie scattering theory, where the Rayleigh approximation (Equation 4) is no more valid, P + (t) is found to increase with D according to a more complicated power law formulation. Following the scattering model of Gouhier and Donnadieu (2008), and for blocks larger than 9 cm, the timeintegrated mass M (kg) becomes: where h is the vertical length of the given probed volume (in m). As in Equation 5, the first factor can be grouped into a constant C ′ (kg mW −1 m −1 ). While most parameters are known (γ , θ ) or could be roughly estimated (ρ, S), the radar-sensitive mean diameter D can hardly be estimated, especially during an eruption, despite its dominant weight in the relationship owing to the power law. For this reason, our approach aims at calibrating the constants C (and hence C ′ ) against results from other methods in order to obtain absolute radar-derived mass load, as explained later. Most interestingly, the integral factor represents a Radar proxy for the mass of tephra M * depending only on radar power and velocity measurements. A proxy for the total erupted mass of tephra (TEM * ) can be obtained by integrating the radar mass proxy M * over the total duration of the paroxysm and over all range gates capturing ascending tephra above the crater. It is also straight forward to compute an average Mass Eruption Rate proxy (MER * ) from the number of samples n recorded at acquisition intervals dt = 0.23 s in a probed volume above the emission source: Time series of the radar-derived MER proxy can thus be computed at high rate (MER * ), potentially at each acquisition time (i.e., at rate 1/dt using n = 1), to inform in real-time on the eruption intensity evolution including during overcast weather preventing visual observations. MER * thus provides an useful tool to improve the real-time monitoring and forecasting volcanic ash dispersal and fallout during lava fountain paroxysms at Etna. Because the range resolution (150 m) is usually smaller than the lava fountain width, several range gates commonly dominate the echoes amplitude and are used for the aforementioned spatial integration of the mass or MER proxy: 3,135 and 3,285 m for the NSEC paroxysms, 3,885, 4,035, and 4,185 m for the VOR paroxysms in December 2015, and 4,035 and 4,185 m for the Northeast Crater. These are the range gates above the erupting crater, as seen from the sounding geometry (Figures 1B,C).
In the next section, we illustrate the use of the radar mass proxy to infer on the dynamics during an explosive eruption at Etna.

Eruption Dynamics During Etna's Paroxysmal Activity
First, paroxysmal eruptions at Mount Etna present a similar succession of eruptive phases (Bonaccorso et al., , 2013Behncke et al., 2014;. The first phase corresponds to a discrete Strombolian activity lasting hours to several days (Behncke et al., 2014), which is not well captured by the radar at the very beginning of the paroxysm owing to tephra emissions mostly confined inside the crater and the lack of sustained plume above the crater rims.
Secondly, the number and intensity of explosions increase and a transition toward an unsteady lava fountain regime occurs (phase 2 in Figure 3, 15-20 min). This period of increasing intensity might represent the evacuation of the partially degassed conduit magma from the previous eruption as it becomes pushed out of the conduit  and replaced by newer magma richer in gas. As new magma progressively fills up the entire conduit, the flow regime transitions from slug flow to churn flow leading to an unstable lava fountaining (Ulivieri et al., 2013). This unsteady phase can be characterized by a shoulder (first bump) in the radar signals, well observed during strong paroxysms like those on 23 February 2013 and 23 November 2013 for example (Figures 3B,C).
Then, two main types of paroxysms can be distinguished during the following third eruptive phase. The 27 Type-A paroxysms (57.4% of the total) are characterized by a clearly sustained climax phase lasting 44.19 ± 5.30 min in average (Figures 3B,C and Table 1). In type-B paroxysms (42.6% of the total), the climax phase is not always well defined, suggesting a lava fountain regime remaining unstable ( Figure 3A). Over 20 Type-B paroxysms, only 8 (40%) present identifiable sustained phases during 44.25 ± 17.89 min in average (Type-B1, Table 1). The 12 (60%) other paroxysms are characterized by highly variable tephra emission (Type-B2).  Finally, the fourth phase (Figure 3) is characterized by a relatively rapid decrease in the radar signal (between 7 and 70 min) with respect to the eruption duration, with an average of 25.4 min during which the lava fountain stops and is replaced by ash emission not well captured by VOLDORAD 2B. Four long-lasting paroxysms present outlier values of decrease time lasting 126 and 289 min (episodes E20, E40, E41, and E43 in Appendix 1).
Average velocities during the climax phase range between 55 and 200 m s −1 (Appendix 1), with a mean of 125 ± 6 m s −1 . However, ejection velocities can reach much higher velocities for a few seconds, 360 m s −1 at the highest over all the paroxysms (short peaks up to 432 m s −1 ). Ejection velocities exceeding 400 m s −1 had previously been measured at Etna using the same type of radar during the Laghetto eruption in July 2001 (Donnadieu et al., 2005). Maximum velocities measured by radar are generally higher than those estimated from infrared. However, in Figure 3, it is important to notice that the pyroclastic emission during a paroxysmal activity is highly variable as a function of time, and variations are also different among the paroxysms. Indeed, ratios shown in histograms of Figures 5A,B indicate that the climax most frequently releases about 80 to 90% of the TEM (modal class) with an average of 76%. Likewise, the MER * averaged over the whole paroxysm most frequently represents about 43% of the climax MER * (Figure 5 and Table 1). Paroxysms after October 2012 show average mass parameters (i.e., TEM * and MER * in Figures 5A,B), during the climax and during the total duration of the events, about twice the averages between 2011 and 2012. This can be a result of a better beam sampling of the lava fountains after the antenna rotation toward the NSEC in October 2012 (Figure 5A,B). Nevertheless, both mass load parameters are homogeneously distributed over nearly two orders of magnitude, indistinctly before (2011-2012) and after (2013-2015) the antenna rotation, with good correlations (e.g., R 2 of 0.98 and 0.94 in Figures 5A,B).

Plume Top Height and Radar Mass Proxy
Plume top heights are strongly controlled by the MER and cross-winds (Morton et al., 1956;Sparks et al., 1997;Bursik, 2001;Mastin et al., 2009). Taking advantage of the capacity of VOLDORAD 2B to efficiently monitor the MER at high rate in real-time and given the ample variations in mass eruption rate observed during lava fountain paroxysms at Etna (Figure 3), we here investigate the relationship of plume heights and the radar-derived MER * . Figure 6 shows times series of radar mass proxy and observed plume top height evolution over the course of four paroxysmal episodes: those on 12 August 2011, 12 April 2012, and 23 November 2013 of the NSEC, and that on 3 December 2015 of VOR crater. For the 12 August 2011 event, heights were measured from the visible camera (ECV), from satellite imagery and from radar. As expected, during each paroxysm, plume top height variations closely follow the radar mass proxy.
For the 12 August 2011, 12 April 2012, 23 November 2013, and 3 December 2015 events, the start of the sudden increase in activity leading to the climax phases according to the radar data occurs 15-21 min before the tephra plumes reach their first maximum heights (Figures 6B-D). From this, to reach 6.3 km (Figure 6A), 4.3 km (Figure 6B), 7.8 km (Figure 6C), and 12.8 km (Figure 6D) above the vent, the estimated upward velocities of the tephra plumes are calculated to be 4. 97, 4.84, 8.67, and 10.67 m s −1 at the very beginning of the climax phase, respectively.
For the 12 August 2011 paroxysm, the plume maximum heights increase significantly before the ascending phase leading to the climax seen in the radar mass. This can be due to the lack of momentum in the waxing phase of this particular weak emission bringing the plume to its top height mainly by simple buoyant upraise before the weak climax has started. The antenna azimuth before October 2012 might also have led to incomplete sampling of weak paroxysms such as the 12 August 2011 (TEM * of 3.64 × 10 −7 mW m) compared to stronger ones like the 12 April 2012 (TEM * of 1.52 × 10 −6 mW m).

Tephra Mass Load Estimates
Although temporal offsets and their variation as a function of time remain to be explained in detail in terms of phenomenology and environmental factors, the evolution of the plume top height during a paroxysm appears closely related to the radar-derived MER proxy. Plume height is an essential input to VATDMs in order to assess hazards from explosive eruptions. The implementation of this capacity of VOLDORAD 2B to provide a relative MER (i.e., a proxy) in real-time and at high rate, in addition to ejection velocities, would already be a step forward in the monitoring of Etna.
However, absolute MER estimates derived from the MER * are of even greater added value. According to Equation (5), converting the radar mass proxy into an absolute mass in kg requires knowledge of parameters constitutive of constant C. However, particle diameters near the source are mostly unconstrained. A way to calibrate C is to compare radar MER proxies with mass eruption rates (MER in kg s −1 ) from empirical laws based on correlation with plume top heights (H in km), such as in Mastin et al. (2009): However, the latter equation is based on a dataset that is biased by the high proportion of strong eruptions, which hence suffers from a lack of more frequent and smaller ones (Woodhouse et al., 2013). Thus, the scaling law of Mastin et al. (2009) does not appear best suited to tephra plumes associated with a MER < 10 6 kg s −1 , also more sensitive to atmospheric conditions common during fountain-fed tephra plumes of Etna. Therefore, we secondly compared with the model of Degruyter and Bonadonna (2012) using wind velocity profiles across tephra plume heights: where ρ a0 is the reference atmosphere density (in kg m −3 ), g ′ is the reduced gravity at the source vent (in m s −2 ) and N is the buoyancy frequency (equals to 1.065 × 10 −2 s −1 for a standard atmosphere). α is the radial entrainment coefficient set at 0.1 (Degruyter and Bonadonna, 2012). β is the wind entrainment coefficient. We used β = 0.5, a value that diminishes the error associated with downwind plume trajectories (Aubry et al., 2017). Finally, v is the wind velocity across the plume height z (in m): Vertical wind profiles were taken from radio soundings operated at the LICT station in the Northwest of Sicily (http://weather. uwyo.edu/upperair/sounding.html). Figure 7A shows the 1 min-averaged Radar MER proxies calculated 5 min before the plume maximum height measurements from visible and satellite imagery during 19 paroxysms of NSEC. Quantitatively, a systematic approach is used to calibrate the MER proxies with a constant C (Equation 5). We consider all data points falling within the plume height-MER model domain of Degruyter and Bonadonna (2012, Equation 9) limited by the vertical profile wind conditions of 0 m/s and highest winds at Etna during the 23 February 2013 paroxysms. The constant needed to reach the highest percentage (i.e., 90% in Figure 7B) of data fitting between the models and the MER proxies is equal to 8.25 × 10 14 kg mW −1 m −1 . Altogether, the MER proxies as a function of observed plume top heights by visible, satellite and X-band radar imagery are scattered on either side of the Mastin et al. (2009)'s statistical law (Equation 8, Figure 7A). Although there is a moderate correlation (R 2 = 0.58), a best-fit power law H ∝ MER * 1/4 is found with an elegant power coefficient of 1/4 well-fitting with the theory (Morton et al., 1956).
First order TEMs can be calculated from the calibrated MERs from VOLDORAD 2B data at Etna: they range from 1.34 × 10 8 to 7.92 × 10 9 kg with an average of 1.37 × 10 9 kg, while the climax erupted masses span from 9.82 × 10 7 to 6.49 × 10 9 kg with an average value of 1.28 × 10 9 kg. Given the radar wavelength, estimated TEMs mainly concern lapilli and block/bombs in the eruptive column. Behncke et al. (2014) have reported a NSEC growth between 2011 and the end of 2012 of about 19 × 10 6 m 3 (bulk volume) due to the proximal fallout. The sum of radar-derived TEMs during the same period leads to a total eruptive bulk volume of detected pyroclasts of 16.1 × 10 6 m 3 . The total erupted bulk volume of detected pyroclasts over all 2013 paroxysms is equal to 26.4 × 10 6 m 3 . This value is also similar to the contribution of proximal fallout, building the NSEC between 2013 and 2014, estimated to 27.0 ± 0.8 × 10 6 m 3 (De Beni et al., 2015). The mean particle density of 1,300 kg m −3 taken to calculate such bulk volumes characterizes the mixture of light (410 kg m −3 , Andronico et al., 2015) and dense block/bombs (close to 2,700 kg m −3 , Bonadonna and Phillips, 2003), and light scoriaceous lapilli (about 600 kg m −3 ; Bonny, 2012) emitted during the paroxysms.
In the next section, we discuss the uncertainties related to the radar mass proxy calibration and the potential benefits of its implementation in real-time for operational monitoring of volcanic activity.

Uncertainties and Implications on Mass Load Parameters
The tephra plume radar sampling has changed on October 2012 because of the antenna rotation eastwards. This might have led to mass load underestimates from radar retrievals of the 2011-2012 lava fountains generating vertical tephra columns. Also, the beam sampling suits better the NSEC lava fountains than the December 2015 VOR paroxysms because the Voragine Crater is more offset from the beam axis ( Figure 1B). The above sampling issues could be highlighted by the three data points falling on the results of Equation 9 based on a 0 m s −1 wind profile (Figure 7). Those points correspond to the 9 July 2011, 1 April 2012 NSEC paroxysms (open circles in Figure 7) and the 3 December 2015 VOR paroxysm (black square).
However, the MER * of all events show consistent distribution trends within a range of two orders of magnitude, whichever the active crater and/or the eruptive periods (Figures 5, 7A). This suggests that the difference in sounding conditions is not a major source of error at first order in mass load estimates from radar parameters. This strengthens our radar-derived mass-proxy methodology to quantitatively characterize the lava fountain paroxysmal episodes of Mount Etna and the high variability of their intensity. Specific environmental conditions such as strong cross wind away from the beam axis, or highly fluctuating wind strength/direction, or strongly bent fountain emission might represent a more significant source of error, underestimating the MER, and these cases should be considered with caution when radar monitoring data are used in real-time for hazard assessment.  Degruyter and Bonadonna [2012;noted D&B (2012), black dashed lines] in the limits of no wind and highest wind vertical profile measured during the 23 February 2013 paroxysmal episode (dashed lines). Triangles refer to plume top heights measured by satellite, circles to ground-based camera in the visible (ECV) and squares to X-band weather radar (Vulpiani et al., 2016).  The calibration of the radar-derived MER provides, even in the absence of constraints on plume height, mass loading parameters that could be used by the INGV-OE to routinely initialize VATDMs. Assuming the particle size distribution does not vary significantly among events, MER * can also be used to directly compare the relative intensity of an ongoing paroxysm with previous ones. In addition to the currently implemented automatic detection and warning of onset and ending of a paroxysmal episode, and the real-time provision of near-source ejection velocities, VOLDORAD 2B could now be further used to automatically locate the active crater by means of the range gating and maximum echoes, and to estimate MER of tephra in real-time with high time resolution. The time series of released mass and hence the mass eruption rate show high variability during an event (Figures 3, 5, 6). This highlights the need to take into account the variations of eruption source parameters during the lava fountains of Etna, in particular the mass-loading parameters, in order to better assess tephra plume hazards. The fact that the MER proxies follow closely the variation of the plume top heights reflects the control of tephra plume ascent by the dynamics of the lava fountains and eruptive column (Figure 6).
The MER for instance is known to strongly control plume height (Mastin et al., 2009). Yet, the average MER is often obtained a posteriori and considered constant, being usually deduced from the total erupted mass inferred from post-eruption deposit analyses and eruption duration. As shown in Figure 5B, the MERs corresponding to the whole paroxysm durations are underestimating by a factor 2.6 the climax MER, and hence potentially the maximum plume height derived from deposit analyses, whereas 76% of the TEM in average is released during the climax ( Table 1). Figures 5A,B emphasize the high contribution of the climax phase in terms of tephra mass load, still assuming that the particle size distribution remains the same during an event.
Thus, the main eruption source parameters are available to operationally initialize dispersion models and constantly reevaluate their input parameters. In fact, by not taking into account cross wind considerations, Mastin et al. (2009) results are supposed to underestimate the MERs for a given plume top height (Degruyter and Bonadonna, 2012). The systematic procedure used to infer a calibration constant of 8.25 × 10 14 kg mW −1 m −1 highlights the spread of our data. Shifting the calibration constant value by a factor of 2 would leave only 80% of the data points inside the Degruyter and Bonadonna (2012) model bounds (Figure 7B). In the case of only 50% matching, the constant varies by a factor of 4-6 toward lower and upper bounds, respectively. Hence, the variation of the calibration constant to obtain absolute MERs, in agreement with Equation 8 and 9, is still reasonable compared to the uncertainty of Mastin et al. (2009)'s formulation (a factor of four within a 50% confidence interval). However, there is still no information concerning the grain size distribution inside lava fountains. The coarsest part of the Total Grain Size Distribution released during Etna's paroxysmal episodes, which falls within the first 5 km from the vents, is rarely sampled (Andronico et al., 2014b;Spanu et al., 2016). Hence, the variability in eruption intensity and fragmentation could lead to different values of the calibration constant C, which is also related to the density and size of detected tephra.

Multi-Method Integration
Eruption Source Parameters are essential to initialize VATDMs in order to forecast the impact of tephra plumes and mitigate related risks. The TGSD is a particularly important parameter to estimate and is not provided by our methodology using instead a calibration of the radar mass proxy against other methods. Owing to the scattering theory, electromagnetic methods are mostly sensitive to a given range of particle sizes as a function of their wavelength. Methods such as satellitebased infrared observations, like SEVIRI in the thermal infrared spectral range (Corradini et al., 2016) discriminate fine ash transported in the atmosphere from micron size up to 20 µm (Wen and Rose, 1994). Samples collected in the field are generally upper limited to centimeter-sized lapilli. The X-band weather radar (DPX4; 3 cm-wavelength) in Catania also used to monitor fountain-fed plumes of Etna is mostly sensitive to particles ranging from 25 µm up to lapilli-sized tephra (Marzano et al., 2012). Comparatively, VOLDORAD 2B Doppler radar (23.5 cm wavelength) is mostly sensitive to cm-sized lapilli up to pluridecimetric blocks and bombs. Thus, each technique provides mass load outputs reflecting the mass proportion of the TGSD fraction for which it is the most sensitive. Unsurprisingly, mass estimates should differ among methods, providing for instance TEM underestimates. The mass proportion of each fraction of the TGSD, however, is poorly known. The mass fraction of block-sized particles in the total, despite its presumably very significant proportion (Spanu et al., 2016), is generally not taken into account. This fraction is well captured by VOLDORAD 2B close to the emission source and mass estimates could be largely improved by integrating its measurements. Behncke et al. (2014) estimated the proportion of distal pyroclastic emission from 2011 to 2012 paroxysms of about 3 × 10 6 m 3 . Such a value corresponds to 14% of the total pyroclastic emission building the NSEC (i.e., the paroxysmal proportion being equal to 19 × 10 6 m 3 ). Accordingly, it means that VOLDORAD 2B, by its detection of the lava fountain, is able to detect the maximum total erupted mass from Etna's paroxysmal episodes. Figure 9 shows the reasonable good agreement between the calibrated TEMs from our radar data and the TEMs retrieved from X-band weather radar, ground-based and satellite-based infrared imagery, and from post-eruption deposit analyses.
Data points scatter across the isomass baseline with 71% of estimates by remote sensing methods within a factor of 3 of our calibrated radar TEMs. Data points falling above the isomass baseline correspond to less material detected by VOLDORAD 2B. These include mainly the VOR Crater paroxysms (3-5 December 2013) and, to a lower degree, the 12 January 2011 NSEC paroxysm. However, those paroxysms, as suggested before, could have been less well sampled, and hence their TEM underestimated by VOLDORAD 2B in Figure 9, owing to activity location offset. This means, therefore, that VOLDORAD 2B TEMs are supposed to be always larger than those found by X-band radar, ground-based and satellite-based infrared data.
In addition, TEMs retrieved from post-eruption deposits by Andronico et al. (2014aAndronico et al. ( , 2015 for the 23 November 2013 and the 12 January 2011 are 2 to 10 times weaker than the remote sensing ones. This is probably due to the lack FIGURE 9 | Total Erupted tephra Mass (TEM) at Etna as a function of the calibrated mass proxy M* from the VOLDORAD 2B radar parameters. The isomass baseline appears as a dashed line. Purple triangles correspond to the mass obtained from X-band radar during four paroxysms of the VOR Crater in December 2015 (Vulpiani et al., 2016), X-band weather radar and satellite during the 23 November 2013 paroxysm (Corradini et al., 2016), and from an integration of field, dispersion model, ground-based and satellite data during the 23 February 2013 event (Poret et al., in press). Orange cross corresponds to the TEM retrieved from ground-based infrared imagery on the 12 January 2011 . Blue dots show the results from post-eruption deposit analyses the 12 January 2011 and the 23 November (Andronico et al., 2014a(Andronico et al., , 2015. of tephra sampling in the first 5-6 km from the NSEC. In fact, Spanu et al. (2016) have shown that a lack of sampling inside the first km from the Southeast Crater, after the 24 November 2006 paroxysm, can lead to a loss of 30% of the TEM. Moreover, Andronico et al. (2014aAndronico et al. ( , 2015 do not consider the deposits on the pyroclastic cone that were instead evaluated by Behncke et al. (2014). Mainly for this reason, we suggest that the total mass derived from deposits should rather be called Plume Erupted Mass instead of Total Erupted Mass of tephra, and this, in the case of a paroxysmal activity involving a fountain-fed plume. Given the particle size overlap among methods, the total grain size distribution could be determined through a multi-frequency combination of remote sensing methods and field sampling, and used in dispersal models (Poret et al., in press). Future efforts should aim at this objective. Indeed, a comparison between radar-inferred TEMs and those obtained by post-eruption deposits could be useful to investigate the issue of partitioning of proximal fallout recorded by VOLDORAD 2B relative to the distal ash mass fraction sampled up to 400 km away from Etna (Andronico et al., 2015).
Finally, the methodology of the radar mass proxy could be transportable to other radars used for the monitoring of other volcanoes. In particular, several scanning Weather Doppler radars are able to measure the above-vent radial velocities and echo power, in addition to detect the whole eruptive column and their internal properties (example of the VARR of Marzano et al., 2006). Weather-radar estimates at the source could be improved thanks to our methodology being independent of the detected particle diameters. Moreover, in terms of multi-method integration, our estimates of near-source ejection velocities from the VOLDORAD data base (Donnadieu et al., 2015) could be used to refine as well the DPX4-inferred MER estimates of the Voragine paroxysms, as suggested by Vulpiani et al. (2016).

CONCLUSIONS
Forty-seven out of the 49 paroxysmal episodes of fountainfed tephra plumes produced by Etna between January 2011 and December 2015 were analyzed using the high rate data of the 23.5-cm wavelength Doppler radar (VOLDORAD 2B) monitoring the explosive activity of the summit craters. A methodology has been developed to compute a radar mass proxy, and hence a Mass Eruption Rate proxy. In addition to the estimation of near-source ejection velocities with a high time resolution, the radar mass proxies allow to study the dynamics of Etna's paroxysms. Although there is limitation of the full sampling of lava fountains in 2011-2012 and 2015 because of the detection angle of the radar beam, each derived mass parameter during the climax phases and the total duration of the paroxysms seem correlated, and this, no matter the detection limits. Paroxysmal episodes of Etna present highly variable volcanic emission as a function of time but the tephra mass released during the climax phases most commonly represent more than 70% of the total erupted mass.
By calibrating the radar MER proxy with models relating MER to plume height, TEMs and MERs are found to correlate with TEM inferred from independent remote sensors. Eventually, the developed mass proxy methodology allows the real-time assessment of eruption source parameters at Etna, but also at any volcano monitored by Doppler radars. These eruption source parameters include vent location, event duration, nearsource ejection velocities, MER evaluation and expected plume top height at first order. This could be integrated into the 24/7 procedure operational during volcanic crises at Mount Etna. Given the lack of information on the total grain size distribution, synergetic efforts should now aim to combine sensors working at different wavelength (radars, ground-based, and satellite imagery) with field deposits analyses to refine the MER and complete TEM during the next paroxysmal activity at Etna, as well as to investigate the partitioning between proximal and distal tephra.

AUTHOR CONTRIBUTIONS
VF-L processed the Doppler radar data. VF-L and FD interpreted the data and wrote the manuscript. SS and MP processed the plume top heights data retrieved by visible imagery. AP developed the original physical model behind the mass proxy methodology. SS and MC provided INGV-OE monitoring data and helped with the writing process of the manuscript. PF, CH, YG, and MP were in charge of radar data acquisition and formatting.