Tracking Formation of a Lava Lake From Ground and Space: Masaya Volcano (Nicaragua), 2014–2017

2.1 Masaya Volcano

Masaya is a tholeiitic basaltic shield volcanic complex (Walker et al., 1993) in Central Nicaragua (Figure 1a), world-renown for its nearly persistent open-vent degassing activity (Stix, 2007) and for having produced some of the few examples of basaltic plinian eruptions (Kutterolf et al., 2007; Pérez et al., 2009; van Wyk de Vries, 1993; Williams, 1983). The Masaya volcano in a strict sense is the youngest shield edifice of the complex and sits on an older sequence of nested calderas and craters. Masaya is cut on its summit by a NW-SE elongated caldera that is 11 km long by 6 km wide (Figure 1b) and was formed by at least three major (plinian) basaltic eruptions between ∼6 and 1.8 ka ago (Bice, 1985; Kutterolf et al., 2007, 2008; Pérez et al., 2009; Pérez & Freundt, 2006; van Wyk de Vries, 1993; Williams, 1983). The post 1.8 ka activity has been concentrated in the caldera center along a semicircular set of vents/cones, today topped by four pit craters, among which is the currently active Santiago crater (Rymer et al., 1998; Figure 1c).

The Santiago pit crater (Figures 1d and 1e) formed in a sequence of collapse/explosion events in 1858–1859 (McBirney, 1956; Rymer et al., 1998) and has since remained the main site of activity at Masaya. During its short-lived history, Santiago has alternated between (i) periods of intense degassing associated with opening of active vents on the crater floor eventually culminating with lava lake formation (such as in 1965–1969, 1972–1979, 1989, and 1993–1994) and (ii) phases of reduced activity. Each “degassing crisis” has typically lasted years to decades (the most recent started in 1993) and has typically fluctuated between vent-opening phases (mostly due to unroofing of small chambers beneath the crater floor) and vent-closing phases (due to rockfalls inside the crater; Rymer et al., 1998). Occasionally, the latter crater wall collapse events have been followed by vent-clearing ash explosions, such as in 2001 (Duffell et al., 2003) and in April–May 2012 (Global Volcanism Program, 2016; Pearson et al., 2012). After a period of reduced volcanic activity between late 2012 and late 2015, incandescence was again reported on the crater floor by INETER on 11 December 2015 (Global Volcanism Program, 2016). Between mid-December 2015 and March 2016, two to three incandescent vents (Figure 1f) became visible, which progressively widened to merge into a single, vigorously active lava lake (Figure 1g), perhaps due to collapsing of the crater floor. Intense seismicity and SO₂ degassing have been associated with formation of the lava lake (Global Volcanism Program, 2016). The lava lake is still active at the time of writing.

2.2 Gas Measurements

We field deployed a Multi-GAS on Masaya in March 2014, as part of the global network of the Deep Earth CArbon DEgassing (DECADE; https://deepcarboncycle.org/home-decade/) project funded by the Deep Carbon Observatory (DCO; https://deepcarbon.net/). This fully autonomous gas sensing unit was installed on the outer rim of the Santiago pit crater (coordinates: 11°59′14.410″N, −86°10′3.734″W, see Figure 1c). The Multi-GAS is still operating by November 2017, but 3 years of results (March 2014 to March 2017) are discussed here. The instrument measured in situ the mixing ratio of CO₂, SO₂, and H₂S in the volcanic gas plume, using the established combination of Near Dispersive Infra Red (NDIR) spectrometers and electrochemical sensors described in previous work (Aiuppa et al., 2014, 2017). Ambient pressure, temperature, and relative humidity were also comeasured, from which the in-plume H₂O mixing ratio were obtained using the Arden Buck equation (Buck, 1981). Data acquisition was controlled by a Moxa embedded computer (model 7112plus), and the whole system was powered by two batteries and a 120 W solar panel. The Multi-GAS unit was programed to acquire data (at 0.1 Hz) during four daily measurement cycles, each 30 min long. The acquired data were stored in an onboard data logger and telemetered to the INETER base station in Managua using a radio link. Concentration data were then postprocessed using the Ratiocalc software (Tamburello, 2015) to obtain time series of the volcanic gas CO₂/SO₂ and H₂O/CO₂ molar ratios. The data set is accessible at (https://doi.org/10.1594/IEDA/100651), which lists (i) raw (mixing ratio) data and (ii) time-averaged CO₂/SO₂ and H₂O/CO₂ ratios for 30 min-long Multi-GAS acquisition windows. No ratio was derived for acquisition windows in which the maximum SO₂ concentration was <3 ppmv threshold and/or low correlation coefficients (R² ≤ 0.6) between gas species were observed. Based on laboratory tests with standard gases, errors in derived ratios are typically ≤15% and ≤30% for CO₂/SO₂ and H₂O/CO₂ ratios, respectively. The sensors’ calibration was also regularly checked in the laboratory using calibration gases at the beginning and end of each measurement subinterval P1–P5 (see below), and no substantial drift or change was observed.

The volcanic CO₂ flux is calculated from combination of coacquired CO₂/SO₂ ratios and SO₂ fluxes. SO₂ flux is obtained for the period March 2014 to February 2017, from evaluation of data of one scanning-DOAS station from the NOVAC network, located ∼1.5 km WSW of Santiago crater at an altitude of 387 m asl (see Figure 1c). The instrument points to the Santiago crater and scans the plume over a flat surface from horizon to horizon, at steps of 3.6°, recording at each step spectra of solar scattered radiation in the 280–450 nm wavelength range (effective range between 300 and 360 nm due to atmospheric and optical filters). The system is a standard NOVAC-Mark I instrument, which specifications are described by Galle et al. (2010). Each scan is completed in about 5 min during daylight. By integrating the column densities of SO₂, derived by the DOAS method (Platt & Stutz, 2008), over a cross section of the plume, and multiplying by the transport speed, the flux of SO₂ is derived. For this evaluation, we used wind speed from the ECMWF ERA-Interim database (Dee et al., 2011) with a spatial resolution of 0.125° × 0.125° and a temporal resolution of 6 h, which was then interpolated to the coordinates of the crater and time of each measurement. Plume direction was obtained by observing the distribution of column densities on each scan, and assuming that the plume centre of mass was drifted at the summit level. Only measurements that were elevated (to ensure complete plume coverage and availability of a background spectrum for cancellation of instrumental and atmospheric effects) were used for further analysis. A total of 28,754 flux measurements were obtained (for details see https://doi.org/10.1594/IEDA/100672), each with an estimated uncertainty of <26% (after Galle et al., 2010). This uncertainty mostly comes from uncertainty in wind speed data, since radiative transfer effects are likely less important due to close proximity to an ash-free plume. Available SO₂ flux data were averaged to obtain mean values for each valid Multi-GAS temporal window. The cumulative error in derived CO₂ fluxes is estimated at 40% from propagation of Multi-GAS and NOVAC uncertainties.

2.3 Space-Based Thermal Data

Thermal data acquired by MODIS sensors were analyzed using MIROVA, an automated, near-real time volcanic hot spot detection system, developed at the University of Turin (Coppola et al., 2016b; www.mirovaweb.it). The MIROVA system uses the MIR radiance measured by the two MODIS sensors (carried on Terra and Aqua NASA's satellites, respectively) that scan the entire earth surface approximately four times per day (2 nighttime and 2 daytime) with a nominal ground resolution of 1 km. These features, together with the presence of a low-gain middle infrared (MIR) channel (MODIS band 21), make MODIS particularly suitable for near-real time monitoring of worldwide volcanic activity (Coppola et al., 2013, 2016b; Rothery et al., 2005; Wright et al., 2002, 2004, 2008).

The hot spot detection algorithm consists of contextual spectral and spatial principles specifically designed to efficiently detect small-scale to large-scale thermal anomalies (from <1 MW to >40 GW) while maintaining false detections low (<1% according to Coppola et al., 2016b). This procedure allows to track a large variety of volcanic activity, including lava flows (Coppola et al., 2017a, 2017b) and domes (Coppola et al., 2015; Werner et al., 2017), as well as strombolian activity (Coppola et al., 2014). Nonetheless, MIROVA has shown its efficiency in detecting previously unknown effusive activity at remote volcanoes (i.e., Nevados del Chillan; Coppola et al., 2016c), in tracking the development of intense fumarolic activity at Santa Ana volcano (Laiolo et al., 2017) and in detecting the recent rebirth of Nyamulagira lava lake (Coppola et al., 2016a).

The Wooster et al. (2003) formulation is used to retrieve the volcanic radiant power (VRP; W) from MODIS data, starting from hot spot pixels detected by MIROVA:

(1)

where L_MIR and L_MIRbk are the MIR radiances (W m⁻² sr⁻¹ µm⁻¹) characterizing the single hot spot pixel and the background. The coefficient 1.89 × 10⁷ (m² sr µm) results from best fit regression analysis between above background MIR radiance and radiant power (Wooster et al., 2003) and allows estimations of VRP (±30%) of hot surfaces having temperatures ranging from 600 to 1,500 K.

Visual inspection of acquired images is used to discard false alerts, isolated fires, thermal data contaminated by clouds, or images affected by poor viewing geometry conditions. Between 2015 and 2017, MIROVA identified 282 hotspots over Masaya volcano. Of these, 8 hotspots were due to false detections, and 31 were associated with wild fires. All of these nonvolcanic sources were easily identified because they occurred at more than 5 km from the volcano summit. Of the remaining 243 thermal anomalies associated with volcanic activity, 88 images were discarded because of the presence of clouds and/or because the high satellite zenith angle (i.e., >45°) impeded a clear view of the Masaya summit craters. The supervised VRP measurements collected during the analyzed period are listed in the supporting information Table S1.

4.1 A CO₂-Rich Gas

Since onset of the 1993 to present degassing unrest, and prior to our study, the H₂O-CO₂-SO₂ signature of the Masaya volcanic gas plume has been determined in six occasions during field surveys in 1998, 1999, 2000, 2001, 2006, and 2009 (see Martin et al., 2010 for a review). These previous studies indicated plume CO₂/SO₂ ratios of 1.5–3.5 (mean: 2.5 ± 0.7) and a hydrous gas signature (H₂O/CO₂ ratios of 27.9 ± 11.2) (Figure 4a).

Our systematic gas observations demonstrate a systematically more CO₂-rich gas vented by Masaya in 2014–2017 (Figure 4a). The time-averaged (all data, 2014–2017 period) CO₂/SO₂ ratio is 6.2 ± 3.5 (range 2.2–41), well beyond the 1998–2009 range. Our derived H₂O/CO₂ ratios are also lower (2014–2017 average: 10.3 ± 8.8). An especially CO₂-rich composition (Figures 2 and 4a) is observed in mid to late-November 2015 (our subinterval P3), a few weeks prior to “rebirth” of the Masaya lava lake in December 2015. In the same subinterval P3, the CO₂ flux (Figure 3c) is 2–3 times the average 2014–2017 levels (26.2–37.3 kg/s), and 7.5 times the Martin et al. (2010) estimate of ∼10.7 kg/s (in March 2009). After P3, gas parameters (composition and fluxes) are back to P1–P2 levels. In summary, our results point to a perturbed degassing regime at Masaya in 2014–2017 and suggest unusually elevated CO₂ release prior to appearance of the lava lake in late 2015.

4.2 Measured Versus Modeled Gas Compositions

The time-evolving volcanic gas composition in 2014–2017 and especially the elevated CO₂ degassing phase P3 are quantitatively interpreted from results of model runs performed with a volatile saturation code (Figure 4). As in previous work (Aiuppa et al., 2007, 2010, 2016; de Moor et al., 2016), the model of Moretti et al. (2003) and Moretti and Papale (2004) is used to simulate the pressure-dependent compositional evolution (in the C-O-H-S system) of the Masaya magmatic gas phase formed upon magma decompression in the 400–0.1 MPa interval (Figure 4b).

The model runs are initialized at conditions representative of the Masaya magmatic system (Table 2). Temperature is set constant throughout the runs at 1,116°C (from de Moor et al., 2013, calculated from the olivine-liquid geothermometer of Putirka (2008)). Oxygen fugacity of Masaya melt is fixed at 1.7 Log units above the QFM (ΔQFM = +1.7, whereby QFM is the quartz-fayalite-magnetite buffer of Frost (1991)), based on measured iron speciation (de Moor et al., 2013). Melt composition and total volatile contents (Table 2) are inferred from compositions of Masaya melt inclusions (MIs; Atlas & Dixon, 2006; de Moor et al., 2013; Sadofsky et al., 2008; Wehrmann et al., 2011; Zurek, 2016) and glasses from laboratory degassing experiments performed at Masaya-like conditions (Lesne et al., 2011).

Melt inclusion results suggest that Masaya volcano is fed by relatively water-poor (1.5–1.6 wt %) magma, similar to other primitive magmas in Central Nicaragua (e.g., Granada and Nejapa volcanoes; Wehrmann et al., 2011). Available melt inclusion data for Masaya show therefore no evidence of the water-rich (H₂O = 3.0–6.1 wt %) magma component seen at other Nicaraguan volcanoes, including Cerro Negro, Telica, San Cristóbal, and Cosigüina (Longpré et al., 2014; Portnyagin et al., 2014; Robidoux et al., 2017a, 2017b; Roggensack et al., 1997). Note that although the degassing experiments of Lesne et al. (2011) were initially prepared and run at 1.5–1.7 wt % H₂O, final total H₂O contents in experimental charges (fluid + glass) were >2.6 wt % because of water production due to hydrogen exchange and reduction of ferric iron (Lesne et al., 2011; Table 2). The experimental results of Lesne et al. (2011) are unique for Masaya, but their higher (than MIs) H₂O contents are likely to have determined vapor-melt partitioning unrepresentative of natural conditions (see below). The Lesne's et al. (2011) experiments will be thus considered in the following for comparative purposes only, but not as stringent constraints for our interpretation.

Melt inclusion information also converges to indicate a preeruptive S content of ∼500 mg kg⁻¹ (Table 2). This implies that the MAS.1.A experiment of Lesne et al. (2011) (our run#2) is likely the closest to natural conditions (apart from the aforesaid issues on the total water content), while the MAS.1.B experiment of Lesne et al. (2011) (our run#3) should only be viewed as a hypothetic S-rich example. As commonly observed at arc volcanoes (Wallace, 2005; Wallace et al., 2015), the preeruptive melt CO₂ content is the least well constrained. Wehrmann et al. (2011) reported one single inclusion with a measured CO₂ content of 369 mg kg⁻¹, and this is used as the basis for our run#1 simulation. In run#4, we consider a CO₂-richer magmatic content, based on the maximum melt inclusion CO₂ content (∼6,000 mg kg⁻¹) of Atlas and Dixon (2006). Lesne et al. (2011) assumed a similarly high CO₂ content (∼7,000 mg kg⁻¹; adopted in our runs#2 and #3).

Model results (Figure 4) show that, at any explored condition, the magmatic gas in equilibrium with the Masaya melt is rich in CO₂ at high pressure. As pressure decreases, the gas composition becomes H₂O-dominated and richer in S in all model runs (Figure 4). The different model trends obtained in the four runs reflect the distinct initial volatile contents used, with the most hydrous gas compositions obtained in the CO₂-poor run#1 and the least hydrous in the S-rich run#3. The four model lines together identify the compositional field of Masaya magmatic gases (Figure 4a), which also encompass the inferred compositional range of fluids formed in the experiments of Lesne et al. (2011).

The model runs also confirm the strong pressure dependence of the CO₂/S ratio (Figure 4b), as seen in other mafic systems (Aiuppa et al., 2007, 2010, 2017; de Moor et al., 2016). The model-derived CO₂/S versus pressure trends are less steep, and shifted toward CO₂-enriched compositions, relative to the experimental trends obtained by Lesne et al. (2011). This discrepancy, particularly evident at high pressure, reflects the 80% higher total water contents in the experiments of Lesne et al. (2011) (see above), which may have acted to determine a gas phase enriched in H₂O and S (at the expense of CO₂). We additionally argue that the CO₂/S ratios in the experiments of Lesne et al. (2011), being relatively constant over a wide P-range (Figure 4b), are inconsistent with a sulfur dissolution behavior simply determined by gas-melt partitioning. Rather, they suggest that the sulfur saturation content at sulfide saturation (SCSS) was likely reached during the experiments, leading to formation of solid/liquid iron sulfides (Mavrogenes & O'Neill, 1999; Moretti & Baker, 2008). Sulfide saturation, while not reported by Lesne et al. (2011), is well consistent with the hydrous nature of their experimental melts (Fortin et al., 2015; Moune et al., 2009). H₂O-rich conditions determine a high Fe²⁺/Fe_tot ratio (>0.7; Lesne et al., 2011 and our modeling) even at oxidation states of NNO + 1 or more (Moretti & Baker, 2008). In contrast, no sulfide separation is predicted by our model runs (conducted at the less hydrous conditions of Masaya natural melts), justifying the disagreement between model (this study) and experimental (Lesne et al., 2011) results at high pressures (Figure 4b). Model and experimental trends converge at low pressure (Figure 4b).

Comparison between models and observations (Figure 4a) suggests that the CO₂-rich gas detected in subinterval P3, prior to emergence of the lava lake (Figure 2e), cannot originate from a shallow source, e.g., from near-surface (0.1 MPa pressure) gas-melt separation. The four most CO₂-enriched P3 gases in Figure 4a, in particular, have CO₂/SO₂ ratios of 24.1–40.8. Depending on the model run used, these correspond to model magmatic gas compositions at respectively ∼9–25 and ∼15–35 MPa pressure (Figure 4b). Note that even higher equilibrium pressures (∼100–200 MPa; Figure 4b) would be obtained using the MAS.1.A experimental fluid line of Lesne et al. (2011). In contrast, the CO₂-poorer gases observed in subintervals P1–P2–P4–P5 mostly fall within the low-pressure (0.1–10 MPa; Figure 4a) domain of the Masaya magmatic gas compositional field, indicating more shallow origin. A cluster of P1–P2–P4–P5 gases, finally, plots outside the Masaya magmatic gas field, and is displaced toward the H₂O corner. Most of the 1998–2009 Masaya gas data are similarly more H₂O rich than the predicted (model) magmatic gas compositions. We ascribe these hydrous gas compositions to entrainment of nonmagmatic H₂O in the plume. Interestingly, these H₂O-rich compositions dominate in subintervals P1, P2, and P5 (mostly referring to the Masaya wet season), while they are typically missing during subintervals P3 and P4 (Figure 2f) (encompassing a typical Masaya dry season). We argue that sampling of recirculated meteoric water, perhaps emitted by low-temperature meteoric-H₂O-dominated hydrothermal fumaroles along the inner crater's walls, may justify the nonmagmatic H₂O-rich gas compositions in subintervals P1, P2, and P5.

4.3 A Deep Trigger?

The above model results imply that the CO₂-rich signature of P3 gas is inconsistent with a shallow magma source. To preserve equilibrium compositions corresponding to ∼9–35 MPa pressure, the source magmatic gas phase must have last equilibrated with (and separated from) the melt at equivalent depths of 0.36–1.4 km (calculated using a magma density of 2,600 kg m⁻³ for Masaya; Stix, 2007). This depth interval corresponds to the roots of the shallow magma reservoir identified by Rymer et al. (1998) and Williams-Jones et al. (2003), based on modeling of gravity data. According to these authors, a shallow magmatic body (<∼1 km) is located directly beneath the currently active Nindiri cone (Figure 5), and the temporal changes in its degree of vesiculation (and thus density) govern the observed gravity variations at Masaya. Vesicularity of the shallow resident magma would, in turn, be modulated by temporal variations in gas bubble supply from deeper convecting magma (Stix, 2007), with a larger bubble influx resulting in periods of decreased gravity and elevated degassing at the surface (e.g., a degassing crisis; Delmelle et al., 1999; Williams-Jones et al., 2003).

Our observations here suggest an increasing supply of CO₂-rich gas bubbles to the shallow magmatic body (Rymer et al., 1998) in November 2015 (Figures 3c and 5b). We propose this (factor ∼3) higher than usual gas bubble supply is sourced by degassing of faster convecting magma, at minimum equivalent depths of 0.36–1.4 km (Figure 5b). We infer below the increase in magma convection rate that is required to justify elevated degassing during P3.

The rate at which magma is convectively rising and degassing (magma degassing rate; Table 1) is initially calculated from the measured SO₂ flux, as in Williams-Jones et al. (2003) and Stix (2007). In these calculations, we use a 2,600 kg m⁻³ magma density (Stix, 2007) and assume complete degassing of a magma with 500 mg kg⁻¹ initial S content (Table 2). We neglect crystal content owing to the poorly porphyric (≪10%) nature of Masaya magmas (Zurek, 2016). Using these numbers, the mean SO₂ flux of 9.5 ± 3.2 kg s⁻¹ converts into a magma degassing rate for subinterval P3 of 3.6 ± 1.2 m³ s⁻¹ (Table 1). This calculation implies a ∼20% increase in magma transport rate, relative to P1–P2–P4–P5 (mean 3.1 m³ s⁻¹; range 3.0–3.3, Table 1) is required to account for the elevated SO₂ degassing in P3. Integrated over the 26 days of duration, we infer that ∼8 Mm³ must have completely degassed during P3 to account for surface SO₂ emissions, or ∼1.2 Mm³ in excess to what would have degassed in the same time interval at the background rate of 3.1 m³ s⁻¹ (Table 1).

We also argue that, during P3, the SO₂ flux increases by only ∼17% (relative to P1–P2–P4–P5), while the CO₂ flux increases by a factor 157% (Table 1). As such, the CO₂ flux would in principle be more appropriate to confine the magma degassing rate increase. However, this operation is complicated by the total CO₂ content in Masaya magmas being poorly constrained (Table 2). Assuming a total (parental melt) CO₂ content of 6,000 mg kg⁻¹, the P3 CO₂ flux would imply magma degassing rate of 5.2 ± 2.6 m³ s⁻¹, higher than the 3.6 ± 1.2 m³ s⁻¹ obtained above from SO₂. This mismatch is reconciled by assuming (i) a higher parental melt CO₂ content (8,600 mg kg⁻¹ of CO₂ are required to obtain a 3.6 m³ s⁻¹ magma degassing rate) or (ii) that there is deeply degassing magma that contributes lots of CO₂ but little SO₂ to the shallow convecting part of the plumbing system, or (iii) by considering that only a fraction (∼70%) of total S (500 mg kg⁻¹) is actually degassed, the remaining fraction remaining dissolved in the melt. The latter is well consistent with the relatively deep (∼9–35 MPa pressure) gas source in subinterval P3, and considering that the model melts still contain 295–462 mg kg⁻¹ S in dissolved form at 10 MPa (model runs#1, #2, and #4). Although we cannot entire resolve between the three hypotheses above, owing to the lack of more precise knowledge on parental melt CO₂ content, we still consider our “deep degassing” hypothesis (a combination of (ii) and (iii) above) well motivated. The CO₂ flux-based magma degassing rate of 5.2 m³ s⁻¹ would imply a degassing magma volume (in 26 days) of 11.7 Mm³, and an excess magma degassing magma volume of ∼7.2 Mm³ (Table 1)

In view of the above calculations, we consider likely that degassing of an unusually high (8.2–11.7 Mm³) magma volume, circulating at >∼9–35 MPa pressure (0.36–1.4 km depth) at an average rate of 3.6–5.2 m³ s⁻¹, sustained elevated CO₂-rich gas supply during the 26 days of subinterval P3. Interestingly, our results are consistent with a general inflation of the Masaya edifice since late 2015, interpreted as due to volume change/pressure increase at ∼2.3–3.8 km beneath the surface (Stephens et al., 2017; C. Wauthier, 2017, personal communication).

4.4 Formation of the Lava Lake

An elevated gas bubble supply from depth may impact stability of the shallow magma reservoir that feeds the lava lake in many different ways. According to Stix (2007), the long-lived (1993 to present) degassing unrest at Masaya, occurring with essentially no eruption of magma, can be explained by either one of the two models below (or from a combination of the two). In the conduit convection model (Kazahaya et al., 1994; Stevenson & Blake, 1998), persistent open-vent activity is thought to be driven by density contrast between continuously ascending (and degassing) buoyant gas-rich magma, and degassed (denser) resident magma, being recycled back in the conduit. In the foam accumulation model (Jaupart & Vergniolle, 1989), instead, open-vent activity is controlled by development of a foam layer at the reservoir-conduit discontinuity, in which quiescent leakage from the foam feeds continuous passive gas emissions, while catastrophic collapse of the foam (by reaching a critical thickness) generates explosive magmatic eruptions (not recently recorded at Masaya). In both models, the gas bubble supply is a critical factor. In the conduit convection model, gas bubbles control density of the ascending magma, with an increase in bubble supply leading to greater buoyancy and faster convection. In the foam accumulation model, a larger influx of gas bubbles leads to thickening of the foam layer, increased gas transport from the foam to the conduit (and, hence, enhanced surface degassing), and to increasing magma level in the conduit. A variable gas bubble supply from depth has, for instance, been evoked as the mechanism driving lava lake level fluctuations at Erta Ale volcano in Ethiopia (Vergniolle & Bouche, 2016).

In view of these considerations, it is reasonable to establish a cause-effect relationship between the elevated gas bubble supply (and surface discharge) in subinterval P3, and the ensuing formation of the Masaya lava lake. Although we cannot exclude the concomitant action of additional factors, including gravitational collapse of the unstable crater floor (Harris, 2009; Rymer et al., 1998), we find it unlikely that the elevated degassing pulse in mid to late-November 2015, and the appearance of the lava lake on 11 December 2015, is a mere temporal coincidence. We propose, instead, that a combination of (i) further upward migration of the bubble-rich magma and (ii) rejuvenation of shallow resident magma by deeply sourced gas bubbles, caused the more buoyant column to migrate upward, driving collapse of the crater floor and ultimately manifesting in a vigorous lava lake (Figure 5c). No gas composition information is available for the most vigorous lava lake activity period (PLF; Figure 3a). However, a decreasing CO₂/SO₂ ratio is observed in the last 12 days of subinterval P3 (Figure 2e), and a carbon-poor composition (relative to P3) is observed in subinterval P4 (22 February to 28 March 2016), with the lake still being present and active (Figure 1). Our gas compositional data are therefore consistent with shallowing of the magmatic source (Figure 4).

Enhanced shallow magma convection in late 2015 to early 2016 is clearly supported by the elevated SO₂ fluxes in subinterval PLF (Figures 3a and 3b). The SO₂ flux during PLF averages at 11.4 ± 5.2 kg/s (versus the 8.1 kg/s mean for subintervals P1–P2–P4–P5), and implies degassing of 4.4 ± 2.0 m³ of magma per second (Table 1 and Figures 5 and 6). For the PLF duration of 64 days, this corresponds to convective ascent, complete degassing, and recycling of ∼24 Mm³ of magma. The excess magma volume, e.g., the mass of magma exceeding what would have degassed in 64 days at the background rate of 3.1 m³ s⁻¹, is thus ∼7 Mm³.

In total, we thus infer that ∼7.2 to ∼8.2 Mm³ more magma than during background intervals (P1–P2–P4–P5) circulated (and degassed) underneath Masaya between November 2015 and February 2016 (Table 1). This elevated shallow magma circulation ultimately culminated in formation and growth of the lava lake (Figure 5c).

Further clues on the processes leading to formation of the lake can be obtained by comparing the above SO₂-derived magma degassing rate (Table 1) with the volumetric magma flux inferred from MIROVA-derived thermal data (Figure 6). In fact, while the SO₂ data allows constraint on the rate at which the magma is convectively rising and degassing within the shallow plumbing system ( ), the thermal data can be used to infer the rate at which the magma reaches the uppermost levels of the magma column (i.e., the lava lake), where it radiates and cools before being cycled back (Q_IR).

Following Coppola et al. (2013) this magma flux can be calculated from

(2)

where is the silica content of the erupted lavas (51.5 wt % for Masaya; Walker et al., 1993). This approach and other similar thermal proxies (see Harris & Baloga, 2009) are commonly applied during effusive eruptions in order to estimate the “TimeAveraged Lava Discharge Rate” (TADR) feeding an active lava flow. The typical error associated with satellite-based estimates of magma flux is ±50%, a value comparable with field-based estimates (Coppola et al., 2013; Harris et al., 2007). This uncertainty mostly comes from the implicit assumption of two end-member thermorheological models (defined as hot and cold models by Harris et al. (2007)) that may characterize the spreading and cooling processes of active lava flows. The two models thus provide an upper and lower value of magma flux for any VRP measurement, being representative of large areas/poorly insulated conditions (the hot model – higher radiant density) and small areas/highly insulated conditions (the cold model – lower radiant density). The real volumetric flux may thus vary within these two boundary estimates depending on the local conditions of area and temperature.

In the case of a convecting lava lake, where no magma is discharged from the magmatic system, the thermal proxy (equation 2) can be used to evaluate what would be the equivalent magma discharge rate in the absence of magma recycling (Coppola et al., 2012, 2013). This corresponds to the volumetric magma flux that would need to be supplied to a lava flow in order to radiate an equivalent amount of thermal energy into the atmosphere. In other words, Q_IR provides the rate at which magma is supplied to the surface before being “discharged” (effusive eruptions) or cycled back into the feeding system (lava lakes/open-vent volcanoes).

In this view, we interpret the two magma fluxes ( and Q_IR) as sampling the rate of magma transport at different depths, being the SO₂-derived flux relative to the magma circulation above the exsolution level of SO₂, and the IR-derived flux relative to surficial magma circulation (i.e., few tens of meters at most; Figure 6). The two magma fluxes may be nonsynchronous due, for example, to the different ascent rate of gas and magma within the shallow plumbing system, or due to the delay it takes a degassing-induced collapses in the crater to reveal the magma. However, on a reasonable time scale of few days, the two fluxes can only be equivalent in the limit condition that all magma entering and degassing within the shallow plumbing system ultimately reaches the free air-magma interface.

During subinterval P3, Q_IR is equal to zero (Figures 5b and 6), since no lava lake is exposed at the bottom of the Santiago crater. However, the higher than normal magma degassing rate (see section 9) suggests ascent (and degassing) of a new batch of gas-rich magma in the magmatic system (see above, Figure 5b). We propose that such increased gas supply, and the added magma volume, lead to overpressurization of the shallow magmatic system, causing a gradual rise of the magma column and the subsequence appearance of the lava lake during subinterval PLF (Figure 5c). During this period, the magma degassing rate ( ) is enhanced by ∼0.8 m³ s⁻¹, thus reaching ∼4.4 m³ s⁻¹ (Figure 6). At the same time, Q_IR gradually increases from 0 to ∼0.4 m³ s⁻¹ (about one tenth of the total SO₂-derived flux), suggesting that ≤10% of the degassing magma can finally make its way to the surface in a lava lake (Figure 6). The progressive increase in the magma level continues until March 2016, when Q_IR stabilizes at ∼0.8 m³ s⁻¹ and the magma degassing rate returns back to normal values (i.e., ∼3 m³ s⁻¹) (Figure 6). In these periods (P4 and P5; Figure 5d), the proportion of degassing magma reaching the free surface of the lake increases to ∼20–30%, thus suggesting an improved, but still incomplete efficiency of the plumbing system in transporting the ascending magma up to the free surface. This condition persists to the time of writing, and will likely continue until a further decrease of the magma degassing rate (e.g., in magma feeding into the shallow plumbing system) may cause the lava lake to disappear (as previously seen at Masaya in historical time). On the contrary, the upper limit condition of magma transport efficiency (∼100%, i.e., with all the degassing magma reaching the lava lake) would be hypothetically reached for Q_IR =  ≅ 3.1 m³ s⁻¹, at which the lava lake would be radiating ∼400 MW (equation 2). In view of the rarity of historical eruptions at Masaya, this condition is reached very infrequently. However, a maximum efficiency of transport to the surface could result in further rise of the lava lake level or even eventual production of a lava flow.

4.5 Comparison With Other Lava Lakes

The calculations above suggest an evident unbalance between magma input ( ) and output (Q_IR) at Masaya. We speculate this unbalance may be causing the rarity of overflow or lateral eruptions at this volcano, and perhaps at other systems (e.g., Villarrica and Erebus) that exhibit similarly stable lava lakes. We anticipate that a systematic analysis of versus Q_IR relationships is required to test, and eventually corroborate, our hypothesis, including observations and modeling at those lava lakes characterized by more dynamic eruptive behavior (such as Nyiragongo (Burgi et al., 2014) and Erta Ale (Vergniolle & Bouche, 2016)).

The versus Q_IR unbalance at Masaya requires a very efficient extraction of gas bubbles from the melt, perhaps due to reduced melt viscosity (in the ∼10² Pa s order; Stix, 2007). Gas separation from melt during open-system degassing drives sinking of degassed magma back into the conduit (Shinohara, 2008), thus posing a limit to the surface emplaced magma volume (Q_IR).

The contrasting magma rheologies can also justify the distinct 2015 behaviours of Masaya and Villarrica lava lakes (degassed, crystal-rich basaltic andesites feeding the Villarrica lava lake have inferred viscosity in the ∼10³ Pa s range (Witter et al., 2004), or a factor 10 higher than Masaya). As reported by Aiuppa et al. (2017), a period of elevated supply of deeply sourced CO₂-rich gas bubbles, similar in duration and magnitude to that seen at Masaya (this study), was observed in early 2015 at Villarica. However, the elevated CO₂ degassing phase at Villarrica culminated into a violent (VEI, Volcanic Explosivity Index = 2) paroxysmal explosion, while no eruption was observed at Masaya, only the emergence of a vigorously degassing lava lake. A combination of different initial volatile contents, magma overpressures, and decompression rates, in addition to distinct magma rheologies, may have concurred to explain the contrasting volcanic behaviours. In any case, the “quiet” lava lake emergence at Masaya is clearly indicative of efficient magma degassing.