Periodic outgassing as a result of unsteady convection in Ray Lava Lake, Mount Erebus, Antarctica

TL;DR

This study uses numerical simulations to show that periodic outgassing pulses in Ray Lake, Mount Erebus, can result from internal convective dynamics rather than periodic magma influx, challenging previous assumptions about conduit influence.

Contribution

It demonstrates that internal convection within the lava lake can produce observable periodic outgassing without periodic input from the conduit, offering a new perspective on volcanic surface signals.

Findings

Periodic surface behavior can arise from convective lake dynamics.

Internal convection can mimic conduit-driven periodicity.

Surface signals may not directly reflect conduit processes.

Abstract

Persistently active lava lakes show continuous outgassing and open convection over years to decades. Ray Lake, the lava lake at Mount Erebus, Ross Island, Antarctica, maintains long-term, near steady-state behavior in temperature, heat flux, gas flux, lake level, and composition. This activity is superposed by periodic small pulses of gas and hot magma every 5-18 minutes and disrupted by sporadic Strombolian eruptions. The periodic pulses have been attributed to a variety of potential processes including unstable bidirectional flow in the conduit feeding the lake. In contrast to hypotheses invoking a conduit source for the observed periodicity, we test the hypothesis that the behavior could be the result of dynamics within the lake itself, independent of periodic influx from the conduit. We perform numerical simulations of convection in Ray Lake driven by both constant and periodic inflow of gas-rich magma from the conduit to identify whether the two cases have different observational signatures at the surface. Our simulations show dripping diapirs or pulsing plumes leading to observable surface behavior with periodicities in the range of 5-20 minutes. We conclude that a convective speed faster than the inflow speed can result in periodic behavior without requiring periodicity in conduit dynamics. This finding suggests that the surface behavior of lava lakes might be less indicative of volcanic conduit processes in persistently outgassing volcanoes than previously thought, and that dynamics within the lava lake itself may modify or overprint patterns emerging from the conduit.

Tables

Table 1. Table 1: Model parameters, dimensionless numbers, and their values ; reference values in bold; range of values explored in parentheses. (1) Moussallam et al. ( 2013 ) , (2) Klein ( 2002 ) , (3) Oppenheimer and Kyle ( 2008 ) , (4) Molina et al. ( 2012 ) , (5) Giordano et al. ( 2008 ) , (6) Costa et al. ( 2009 ) .

Parameter	Symbol	Units	Values
Melt Density1	${\rho }^m$	kg/m3	2545
Crystal Density1,2	${\rho }^{\chi }$	kg/m3	2720
Vapor Density3	${\rho }^{\varphi }$	kg/m3	1
Thermal Conductivity4	k	W/m/K	1.53
Bubble Diffusivity	${\kappa }_{\varphi }$	m2/s	10-6
Heat of Crystallization	$L_{\chi }$	J/kg	0, 4.1$\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟓}}$
Heat Capacity4	cP	J/kg/K	1367
Solidus1	Tsol	$\text{celsius}$	884.5
Liquidus1	Tliq	$\text{celsius}$	1030
Crystallinity Exponent	q		1.25
Bubble Radius	$a_0$	cm	2.5, 5, 10
Bubble Hindering	$\mu$		3
Reference Melt Viscosity5	$\eta$	Pa.s	104
Rheology Exponent	n		1, 2, 3
Cohesion	c	Pa	0,${10}^7$
Solid Exponential Weakening6	${\lambda }_s$		4.5
Liquid Exponential Hardening5	${\lambda }_l$		3.5
Critical Crystallinity6	${\chi }_{crit}$	vol	0.6
Crystallinity Step Width6	$w_{\chi }$	vol	0.2
Cooling Time	${\tau }_T$	hr	10,15,20
Outgassing Time	${\tau }_{\varphi }$	s	40, 80, 160, $\infty$
Outgassing Depth	$\delta$	m	0.5, 1, 2
Inflow Temperature	$T_0$	$\text{celsius}$	965, 970, 975
Inflow Vesicularity	${\varphi }_0$	vol	0.1, 0.2, 0.3
Inflow Velocity	$u_{in}$	m/s	0.1, 0.2, 0.3, 0.4, 0.5
Inflow Symmetry			YES, NO
Inflow Period	${\tau }_{u_{in}}$	min	5, 10, 20, $\infty$
Dimensionless Number
Inflow/Convection	$R_{\mathrm{in}}$		5.0 (2.5–12)$\times {10}^{-2}$
Segregation/Convection	$R_{\mathrm{segr}}$		2.2 (0.56–8.9)$\times {10}^{-5}$
Gas Removal/Advection	${\mathrm{Da}}_{\varphi }$		3.9 (3.4–4.4)$\times {10}^{-3}$
Heat Removal/Advection	${\mathrm{Da}}_T$		5.7 (4.3–8.6)$\times {10}^{-6}$
Crystal Reaction/Advection	${\mathrm{Da}}_{\chi }$		3.1 (2.1-5.6)$\times {10}^5$
Gas Advection/Diffusion	${\mathrm{Pe}}_{\varphi }$		8.1 (7.1–9.1)$\times {10}^6$
Heat Advection/Diffusion	${\mathrm{Pe}}_T$		8.1 (7.1–9.1)$\times {10}^7$
Sensible/Latent Heat	$\mathrm{St}$		3.2, $\infty$

Table 2. Table 2: Parameter variations for testing convective regimes . Reference values highlighted in bold.

Parameter	Symbol	Units	Trials	$R_{\mathrm{in}}$ ($\times {10}^{-2}$)	${\mathrm{Da}}_{\varphi }$ ($\times {10}^{-3}$)
Inflow Vesicularity	${\varphi }_0$	Vol	0.1, 0.2, 0.3	9.9, 5.0, 3.3	3.4, 3.9, 4.4
Inflow Velocity	$u_{in}$	m/s	0.1, 0.2, 0.3, 0.4, 0.5	2.5, 5.0, 7.4, 9.9, 12	3.9
Outgassing Time	${\tau }_{\varphi }$	s	40, 50, 80, 160, $\infty$	5.0	7.7, 6.2, 3.9, 1.9, 0

Table 3. Table S1: Mineral densities for the calculation of weighted average crystal density. The crystal volume fractions are an average of measurements from Moussallam et al. ( 2013 ) , and density of mineral phases from Klein ( 2002 ) .

Mineral	Volume Fraction	Density (kg/m3)
Anorthoclase feldspar	0.91	2580
Titanomagnetite	0.03	5180
Olivine	0.02	3800
Clinopyroxene	0.02	3200
Fluorapatite	0.02	3200
Average		2720

Equations

$$\mathrm{Re}=\frac{\rho uL}{\eta },$$

$$\nabla \cdot \mathbf{v}$$

$$-\nabla \cdot 2\bar{\eta }\underline{\mathbf{D}}(\mathbf{v})+\nabla P$$

$$\underline{\mathbf{D}}(\mathbf{v})=\frac{1}{2}(\nabla \mathbf{v}+[\nabla \mathbf{v}{]}^T).$$

$$\bar{\rho }=(1-\varphi )[(1-\chi ){\rho }_l+\chi {\rho }_{\chi }]+\varphi {\rho }_g,$$

$$\chi =\frac{{\varphi }_{\chi }}{{\varphi }_{\chi }+{\varphi }_l},$$

$$\frac{\partial T}{\partial t}=-\mathbf{v}\cdot \nabla T+{\kappa }_T{\nabla }^{2}T+\frac{L_{\chi }}{c_P}{\Gamma }_{\chi }+{\Gamma }_T,$$

$${\Gamma }_T=-\frac{T-T_{atm}}{{\tau }_T}\exp \left(-\frac{z}{\delta }\right),$$

$$\frac{\partial \chi }{\partial t}=-\mathbf{v}\cdot \nabla \chi +{\Gamma }_x,$$

$${\Gamma }_{\chi }=-\frac{\chi -{\chi }^{eq}}{{\tau }_{\chi }},$$

$${\chi }^{eq}={\left(\frac{T-T_{liq}}{T_{sol}-T_{liq}}\right)}^q,$$

$$\frac{\partial \varphi }{\partial t}=-\nabla \cdot \varphi {\mathbf{v}}_g+{\kappa }_{\varphi }{\nabla }^2\varphi +{\Gamma }_{\varphi },$$

$${\mathbf{v}}_g=\mathbf{v}-\frac{2{a_0}^2}{9\bar{\eta }}(1-\varphi {)}^{\mu }\Delta \rho \mathbf{g},$$

$${\Gamma }_{\varphi }=-\frac{\varphi }{{\tau }_{\varphi }}\exp \left(-\frac{z}{\delta }\right),$$

$$\bar{\eta }={\left[{\eta }_s\exp (-{\lambda }_s(1-\chi ))\right]}^X\times {\left[{\eta }_{\ell }\exp ({\lambda }_l\chi )\right]}^{(1-X)},$$

$$X=\frac{1}{2}\left[1+\tanh \left(\frac{\chi -{\chi }_{\mathrm{crit}}}{w_{\chi }}\right)\right],$$

$$x=l_0{x}^{\prime },z=l_0{z}^{\prime },\mathbf{v}=v_0{\mathbf{v}}^{\prime },{\mathbf{v}}^v=v_0{{\mathbf{v}}^v}^{\prime },t=\frac{l_0}{v_0}{t}^{\prime },$$

$$P={\rho }_0{g}_0{l}_0{P}^{\prime },T=T_0{T}^{\prime },\varphi ={\varphi }_0{\varphi }^{\prime },\chi ={\chi }_0{\chi }^{\prime },$$

$$\eta ={\eta }_0{\eta }^{\prime },\rho ={\rho }_0{\rho }^{\prime },\mathbf{g}=g_0\hat{\mathbf{z}},$$

$$v_0=\frac{{\rho }_0{g}_0{l}_0^2}{{\eta }_0},$$

$$\nabla P$$

$$\nabla \cdot \mathbf{v}$$

$$\frac{\partial \varphi }{\partial t}$$

$$\frac{\partial \chi }{\partial t}$$

$$\frac{\partial T}{\partial t}$$

$${\mathrm{R}}_{\mathrm{segr}}=\frac{2a_0^2\Delta \rho g_0}{9{\eta }_0{v}_0}\propto \frac{a_0^2}{l_0^2}.$$

$$\mathrm{D}{\mathrm{a}}_{\varphi }$$

$$\mathrm{D}{\mathrm{a}}_{\chi }$$

$$\mathrm{D}{\mathrm{a}}_{\mathrm{T}}$$

$$\mathrm{P}{\mathrm{e}}_{\varphi }$$

$$\mathrm{P}{\mathrm{e}}_{\mathrm{T}}$$

$$\mathrm{St}=\frac{c_p{T}_0}{L_{\chi }}.$$

$${\mathrm{R}}_{\mathrm{in}}=\frac{u_{in}}{u_{cnvt}}=\frac{u_{in}{\eta }_0}{{\varphi }_0\Delta {\rho }_0{g}_0{l}_0^2},$$

Full text

1 Introduction

Lava lakes provide a rare opportunity for direct observation of the near-surface portion of volcanic plumbing systems. There is a potential, therefore, to use observations at lava lakes to improve our understanding of volcanic activity. Realizing this potential hinges on the ability to distinguish the effects of intra-lake processes from those deeper in the plumbing system on surface observables. The goal of this paper is to advance our ability to interpret surface data from lava lakes by isolating how different convective regimes in the lake are reflected in surface observations.

We focus on the case of Ray Lake on Mount Erebus, Antarctica. Similarly to other persistently active volcanoes (Blackburn et al., 1976; Francis et al., 1993; Huppert and Hallworth, 2007; Palma et al., 2011; Kazahaya et al., 1994; Stevenson and Blake, 1998), Erebus emits orders of magnitude more heat and gas than would be expected based on the volume of lava erupted alone. This imbalance indicates continual recirculation of lava in the conduit with volatile-rich, buoyant magma ascending and volatile-poor, denser magma descending. Compared to other lava lakes, Erebus exhibits an intermediate level of convective vigor. It lacks the organized, well-defined, rigid plates found at Kilauea and Erta Ale (Harris, 2008; Harris et al., 2005; Patrick et al., 2016), but neither does it show a severely disrupted surface indicating chaotic circulation like on Ambrym (Carniel et al., 2003; Lev et al., 2019). In this intermediate convective regime, we hypothesize that both conduit flow and the lava lake convection are reflected in surface observables.

The behavior of Ray Lake is characterized by periodic cycles of increased outgassing with a period of 5-18 min, punctuated by sporadic Strombolian eruptions (Calkins et al., 2008; Oppenheimer and Kyle, 2008; Oppenheimer et al., 2009; Peters et al., 2014a, b; Sweeney et al., 2008). The periodic outgassing is phase locked with variations in surface velocity, temperature, heat flux, and gas composition (Oppenheimer et al., 2009). These observations suggest that the physical transport of magma controls the outgassing, but do not clarify which physical process leads to the periodicity.

Bidirectional flow in conduits is prone to instabilities (Beckett et al., 2014; Huppert and Hallworth, 2007; Kazahaya et al., 1994; Stevenson and Blake, 1998; Suckale et al., 2018), which may explain periodic signals at Mount Erebus (Oppenheimer et al., 2009). By this interpretation, the surface record is assumed to be indicative of processes occurring deeper in the conduit. It implicitly assumes that lake dynamics do not exert strong controls on surface observables. Here, we test this assumption through numerical simulations. We hypothesize that convection within the lava lake itself may lead to the observed periodicity, even in the absence of conduit-driven fluctuations.

We use a model of gas-buoyancy driven convection coupled to evolving vesicularity, temperature, and crystallinity-dependent rheology. First, we impose a constant inflow of bubble-rich magma from the conduit to isolate the effects of lake convection. We test a range of magma inflow, gas loss, and heat loss rates to investigate pertinent convective regimes. Secondly, we consider how lake behavior modulates periodicities imposed from the conduit. We compare model outputs of surface velocity, gas flux and heat flux to observational data. Our results reveal several regimes of lake circulation. We identify their characteristic signatures and map out their controlling parameters. We find that lake convection alone can produce periodic surface behavior similar to observations without requiring periodicity in the conduit.

2 Observational constraints on Ray Lake

Mount Erebus is part of the volcanic group on Ross Island, Antarctica. It supports a persistently active lava lake located at the bottom of the Inner Crater (Fig. 1). Ray Lake has been openly convecting at least since its discovery in 1972, and appears to be in near steady state (Calkins et al., 2008; Dibble et al., 2008; Peters et al., 2014b). Observational studies have employed visual, infrasound, and gas monitoring to characterize the convective and explosive activity (Dibble et al., 2008; Peters et al., 2014a). Averaged over several hours, the system has approximately constant surface temperature, heat flux, and gas flux. Superimposed on this steady mean are cycles of increased outgassing, heat loss, and surface velocities with periods of 5 to 18 minutes (Calkins et al., 2008; Oppenheimer and Kyle, 2008; Oppenheimer et al., 2009; Peters et al., 2014a, b; Sweeney et al., 2008). The steady and periodic behavior are occasionally punctuated by non-periodic Strombolian eruptions in which mild explosions disrupt the surface and partially evacuate the lake (Dibble et al., 2008). Larger explosions have revealed a simple cone-like lake geometry $\sim$40 m across and $\sim$30 m deep, with a 4-10 m wide conduit opening at the base of the lake (Oppenheimer et al., 2009).

Based on thermal imaging (Calkins et al., 2008; Oppenheimer et al., 2009), lake surface temperatures vary from $\sim 275-900\celsius$ with mean temperatures $\sim 525-750\celsius$ (Calkins et al., 2008). Calkins et al. (2008) found periodic increases in maximum and mean temperature of $\sim 120\celsius$ and $\sim 20\celsius$, respectively, every $\sim 5$ min, with thermal maxima lasting $\sim 1-2$ min. Ground-based thermal flux estimates at high lake level (surface area of $\sim$ 1400 m2) are 30 $\pm$ 10 MW (Calkins et al., 2008; Oppenheimer et al., 2009). Satellite observations show higher maximum radiant fluxes up to 100 MW (Wright and Pilger, 2008).

A 5-18 min periodic behavior is evident in the mean and peak surface velocities (Peters et al., 2014a). Mean surface velocities are typically less than $\sim$ 0.1 m/s but can reach up to $\sim$ 0.15 m/s at high activity (Oppenheimer et al., 2009; Peters et al., 2014a). Peak velocities vary between $\sim$ 0.1-0.5 m/s over periods of 5-10 minutes (Oppenheimer et al., 2009) but may exceed 0.8 m/s when the lake is most active (Calkins et al., 2008; Lev et al., 2019). Strombolian eruptions are not correlated with the periodic activity and do not regulate surface velocities beyond the $\sim$5 min in which the lake refills (Peters et al., 2014a).

Gas emission studies show a periodic total flux and composition phase-locked with surface activity (Peters et al., 2014a). Plume measurements of total SO2 flux show a dominant 10 min periodicity with a flux of $\sim 0.7\pm 0.3$ kg/s SO2 (Sweeney et al., 2008). Periodic peaks in surface activity are associated with emission of gas enriched in SO2, H2O, HCl, and HF with respect to CO2 and an increase in CO2/OCS (Oppenheimer and Kyle, 2008; Oppenheimer et al., 2009; Peters et al., 2014a). Calculations based on gas composition and SO2 flux measurements yield an estimate for the total gas flux of $\sim$27.3 kg/s from the volcano summit, although it is uncertain what portion of the total flux derives from passive outgassing apart from lake activity.

3 Methods

We perform numerical simulations of a simplified 2D representation of Ray Lake using a mixed finite-element, finite-difference code (based on Keller et al. (2013). The model takes a phase-averaging approach in which solid, liquid, and gas phases are described as interpenetrating continuum fields at the system scale representing volume-averaged phase interactions at the local scale (Drew, 1983).

In this framework, the magma is an aggregate of liquid melt, solid crystals and vapor bubbles where the collective flow of the aggregate dominates the system-scale dynamics. The phases move collectively with the magma velocity determined by solving the Stokes equations, while bubbles are allowed to segregate vertically based on a hindered-Stokes law. The bubble fraction also diffuses as a result of local-scale fluctuations of bubble motion relative to system-scale flow (Segre et al., 2001; Mucha et al., 2004). We do not allow for segregation of crystals, which is appropriate for small grains with minimal buoyancy contrast to the carrier melt. The relatively high viscosity ($>10^{4}$ Pa s) does not permit significant settling over the timescale of lake convection cycles, even of the 5-10 cm anorthoclase megacrysts present at Erebus (Molina et al., 2012; Moussallam et al., 2013). Additionally, we consider crystals to be in near-equilibrium with a melt of constant composition. Thus, the main effect of crystals in the system is their stiffening effect on the aggregate viscosity, which becomes shear-thinning at high strain rates.

3.1 Model description

We estimate the Reynolds number of convection in Ray Lake using

[TABLE]

where $\rho$ is the liquid density ($\sim$2600 kg/m3), $u$ is the characteristic flow velocity ($\sim$0.2 m/s), $L$ a characteristic length we take as the conduit radius ($\sim$5 m), and $\eta$ the dynamic viscosity ($\sim$104 Pa s), yielding $\mathrm{Re}\approx 0.26$. The low value for $\mathrm{Re}$ justifies neglecting inertial terms.

Mass and momentum are conserved for Stokes flow of an incompressible fluid:

[TABLE]

where $\mathbf{v}$ is the velocity, $\bar{\rho}$ is the density, $\bar{\eta}$ the viscosity and P the pressure in the aggregate volume, g the gravitational acceleration and $\underline{\mathbf{D}}(\mathbf{v})$ the deviatoric strain rate tensor:

[TABLE]

The aggregate density $\bar{\rho}$ is the volume average of phase densities:

[TABLE]

where $\rho_{l}$ is the density of the melt, $\rho_{\chi}$ a weighted average density of the minerals present (Klein, 2002; Moussallam et al., 2013) and $\rho_{g}$ the gas density (Table 1); $\phi$ is the volume fraction of vapor in the aggregate volume, and $\chi$ the volume fraction of crystals relative to the silicate volume:

[TABLE]

where $\phi_{\chi}$ and $\phi_{l}$ are the volume fraction of crystals and melt per unit volume of the aggregate, respectively.

We couple equations 2a and 2b to the conservation of energy,

[TABLE]

where the specific heat capacity, $c_{p}$, and the thermal diffusivity, $\kappa_{T}$, are both assumed constant, $L_{\chi}$ is the latent heat of crystallization, and $\Gamma_{T}$ a parameterized surface cooling rate:

[TABLE]

where $T_{atm}$ is the air temperature above the lake, buffered at the condensation point of water vapor (100 \celsius), $\tau_{T}$ is the characteristic cooling time, $z$ the depth coordinate, and $\delta$ the characteristic length of the cooling boundary layer.

We model the evolution of crystallinity by an advection-reaction equation

[TABLE]

where $\Gamma_{\chi}$ is a volumetric crystallization rate given by:

[TABLE]

with $\tau_{\chi}$ the characteristic time of crystallization adjusted to be rapid compared to advective transport, and $\chi^{eq}$ the equilibrium crystallinity taken as a function of temperature using a power-law fit of the form:

[TABLE]

where $T_{sol}$ is the solidus, $T_{liq}$ the liquidus, and $q$ an exponent near unity. The fitting parameters are determined to approximate the equilibrium crystallinity with $T$ reported by Moussallam et al. (2013).

The vesicularity evolves as

[TABLE]

where $\kappa_{\phi}$ is the diffusivity due to local-scale fluctuations in bubble motion (Mucha et al., 2004; Segre et al., 2001), and $\mathbf{v}_{g}$ is the gas velocity given by a hindered-Stokes segregation law:

[TABLE]

where $a_{0}$ is the average radius of bubbles or bubble clusters, $1\leq\mu\leq 5$ is the hindering exponent Richardson and Zaki (1954), which we set to 3 (Manga, 1996); $\Delta\rho$ is the density contrast between vapor and crystal-melt mixture. In (11), $\Gamma_{\phi}$ is a parameterized outgassing rate,

[TABLE]

in which $\tau_{\phi}$ is the characteristic outgassing time, and $\delta$ the same as the cooling boundary layer depth.

Laboratory experiments show a strong dependence of viscosity on crystal content (see Costa et al. (2009), and refs therein). The particle-stiffening effect was described by Krieger and Dougherty (1959) and numerous follow-up studies as recently reviewed by Mader et al. (2013) in the magmatic context. We assume a viscosity model with a smooth step increase where crystals become closely packed and form a contiguous solid (Fig. S1):

[TABLE]

where $\eta_{l}$ and $\eta_{s}$ are the pure melt and solid viscosities, respectively (Giordano et al., 2008), $\lambda_{l}$ and $\lambda_{s}$ are slopes of $T$-dependent viscosity away from the step, and X is a smooth step function,

[TABLE]

which is centered about a critical crystallinity (Arzi, 1978; Costa, 2005), $\chi_{crit}$, and has a width of $w_{\chi}$ (Costa et al., 2009). We cap the maximum viscosity at 1012 Pa s, which is below the viscosity of solid rock at these conditions but presents a large enough contrast with the lake interior to result in rigid lake walls on the time scale of interest. To limit complexity, we do not include the effects of bubbles on magma viscosity.

Crystal- and bubble-bearing magmas exhibit strain-rate dependent behavior (Saar et al., 2001; Caricchi et al., 2007; Renner et al., 2000; Heymann et al., 2002; Costa et al., 2009; Pistone et al., 2012; Mader et al., 2013; Le Losq et al., 2015). We impose a shear-thinning power-law rheology above $10^{-3}$ s*-1* (Caricchi et al., 2007).

The effects of thermal expansivity and pressure compressibility of phase materials for $\Delta T\approx 50\celsius$, and $\Delta P\approx$ 50 kPa are small compared to the density contrasts between phases, therefore we neglect $P,T$-effects on density. We also find that the bubble segregation velocity, $\mathbf{v}^{v}$, remains small with respect to the magma velocity field, and hence its non-zero divergence may be neglected. By neglecting both gas compressibility and the non-solenoidal flow of multi-phase segregation, we obtain an incompressible flow model.

The imposed parameterized surface outgassing and cooling rates represent the collective effects of diffusive, advective, and radiative heat loss from the lake surface (eqs. 7 & 13). Therefore, we do not resolve small-scale processes contributing to the removal of gas and heat near the surface, which are thought to include bubble coalescence and fragmentation (Blower, 2001). We select the characteristic scales for outgassing and heat removal to be consistent with field observations of gas and heat flux assuming a 1400 m2 surface area (Calkins et al., 2008; Oppenheimer et al., 2009; Wright and Pilger, 2008).

3.2 Dimensional analysis

We identify the characteristic physical scales of the problem by performing a dimensional analysis of the governing equations. We scale variables and parameters by the following dimensional scales:

[TABLE]

where

[TABLE]

is the characteristic speed of a Stokes diapir of radius $l_{0}$, the conduit radius, $g_{0}$ the acceleration due to gravity, $\rho_{0}$ the liquid density, $\eta_{0}$ the liquid viscosity at the inflow temperature, $T_{0}$, $\phi_{0}$ the inflow vesicularity, and $\chi_{0}$ the crystallinity at $T_{0}$. We substitute these scales into the governing equations and drop primes to find the dimensionless form:

[TABLE]

A set of dimensionless numbers emerges from this analysis. We define $\mathrm{R_{segr}}$ as the ratio between the characteristic speeds of bubble ascent and magma convection:

[TABLE]

In our models, $l_{0}$ is much larger than the bubble radius ($2.5$ cm $\leq a_{0}\leq 10$ cm), suggesting that magma convection should prevail over bubble segregation.

We find three Damkoehler numbers which compare the volumetric rates of outgassing, crystallization, and cooling to the rate of magma convection:

[TABLE]

For small values of $\mathrm{Da_{\phi}}$, we expect outgassing to be inefficient and gas to accumulate in the domain, while large values of $\mathrm{Da_{\phi}}$ would result in efficient outgassing. When $\mathrm{Da_{T}}$ is small we expect surface cooling to dominate and a stagnant-lid regime to prevail, whereas for large $\mathrm{Da_{T}}$ we expect convective transport to exceed surface cooling and for the lake to remain open. Here, we use only large $\mathrm{Da_{\chi}}$ consistent with near-equilibrium evolution of crystallinity.

We find two Péclet numbers relating advective and diffusive transport of bubbles and heat:

[TABLE]

Both Péclet numbers are very large ($\mathrm{Pe_{\phi}}>10^{6}$ and $\mathrm{Pe_{T}>10^{7}}$) in our simulations, consistent with the advection-dominated transport of gas and heat expected for Erebus.

We compare the ratio of sensible heat to latent heat via the Stefan number,

[TABLE]

We find $\mathrm{St}\approx 3.2$ for parameters of interest here and hence expect that latent heat may be important with respect to sensible heat.

A final dimensionless parameter, $\mathrm{R_{in}}$, compares the imposed inflow speed at the conduit mouth to the characteristic speed of free convection:

[TABLE]

For high values of $\mathrm{R_{in}}$, the inflow rate exceeds the rate of convective transport, which should result in a pile-up of buoyant material near the boundary. When $\mathrm{R_{in}}$ is small, material is removed from the inlet more quickly than it is fed in, which should result in a dripping instability.

3.3 Numerical model setup

3.3.1 Model discretization

The velocity-pressure solution is obtained by a continuous-Galerkin (CG) finite-element method, with variables discretized on a regular mesh of rectangular elements with linear shape functions for velocity and piece-wise constant ones for pressure (Brenner and Scott, 1994; Keller et al., 2013). To avoid the need for stabilization in a CG-based formulation, we discretize the advection-diffusion-reaction equations for temperature, vesicularity and crystallinity by staggered-grid finite-differences on a collocated grid.

The model uses a second-order accurate, upwind-biased Fromm’s method for grid-based advection (Fromm, 1968). Temperature and crystallinity are advected using the aggregate velocity only, whereas the vesicularity is advected based on the vapor velocity, which combines magma convection and bubble segregation ((12)).

3.3.2 Geometry

A novel component of our model setup is that the specific geometry of the lava lake can evolve with the flow dynamics. As illustrated in Fig. 1, our two-dimensional model domain is rectangular and includes the cold host rock around the lake. The surface extent of the lake and the width of the conduit mouth are initially set at 40 m and 10 m diameter, respectively, consistent with field observations of a cone-shaped lake (Dibble et al., 2008; Oppenheimer et al., 2009). We initiate simulations with an $\sim 1.5$ m thick internal thermal boundary layer that linearly connects the lateral edges of the lake to the conduit. The initial temperature of the lake interior is set to the upwelling magma temperature, T = 970$\celsius$, consistent with the observed $\sim$30% crystallinity in the lake (Moussallam et al., 2013); lake vesicularity is initially set to $\phi$ = 0 everywhere. The walls of the lake are initiated with the same, temperature-dependent material properties as the magma in the lake interior with an initial temperature of T = 100$\celsius$, assuming it is buffered by water vapor condensation. Throughout the simulation, the lake bed can dynamically evolve as a result of diffusive cooling, crystallization and thermal erosion by magma flow.

The side boundaries of the rectangular domain are no-slip (u = w = 0) and insulating (${\partial T}/{\partial\mathbf{n}}$ = 0 normal to boundary). Due to the relative inefficiency of thermal diffusion, the lake walls remain close to their initial T = 100$\celsius$, resulting in a high enough viscosity to keep the walls essentially rigid over a model run. The velocity field in the walls where the viscosity remains at maximum (1012 Pa s) remains fixed at zero and is excluded from calculations for reasons of efficiency. The effect is akin to a no-slip condition along the lake bed.

Except for the conduit itself, the bottom boundary is no-slip (w = u = 0) and isothermal. Within the conduit, in- and outflow of magma is imposed by a vertical velocity of sinusoidal shape: $\sin(2\pi x)$ for asymmetrical flow, and $(\cos(2\pi x)+\cos(4\pi x))/2$ for symmetrical flow. The symmetrical condition is broadly consistent with core-annular flow which has been suggested to dominate bidirectional flow in volcanic conduits (Beckett et al., 2011; Stevenson and Blake, 1998; Suckale et al., 2018). The horizontal velocity remains fixed at u = 0. To avoid forced erosion of the thermal boundary layers along the conduit mouth, the velocity profile is applied on the central 8 m of the 10 m conduit only. We perform simulations with both constant and periodic influx. The latter is captured through sinusoidal variability to the inflow speed with time.

The top boundary is free-slip (${\partial u}/{\partial z}$ = 0, w = 0). Gas and heat are extracted from the top of the lake by the distributed sink terms $\Gamma_{T}$ and $\Gamma_{\phi}$ that decay exponentially away from the surface over the depth $\delta$ (eqs. 7, 13). Additionally, we allow vapor segregation across the lake surface but observe that the hindered-Stokes law alone results in a minor outgassing. Disruption of a chilled surface layer by bulging under local gas accumulation is inhibited by this technique, which biases the simulations towards a stagnant-lid regime. To keep the lake open, we select characteristic cooling times sufficiently long to prevent significant freezing of the surface during the model run time.

4 Lava lake simulation results

4.1 Reference parameters

Observational constraints leave considerable uncertainty in the parameters governing the dynamics of Ray Lake. We choose a set of reference parameter values and associated dimensionless numbers, summarized in Table 1, that represent our best estimate of conditions in Ray Lake.

Figure 2 shows a snapshot of temperature, crystallinity, vesicularity, and density fields from the reference parameter test. The crystallinity directly follows the distribution of temperature and the two are therefore shown together in Fig. 2a. The range of crystallinity remains small within the lake. Hence, it imparts only minor variations to the density field, which is dominated by the effects of vesicularity (Fig. 2b). The density/vesicularity field shows the upwelling, gas-rich magma contrasting with the downwelling, outgassed magma.

Over the duration of the reference simulation, we observe two distinct flow regimes. The first 15 min of the simulation mark a spin-up period in which flow disrupts the initial conditions. The first persistent flow regime is characterized by a dripping instability in which roughly equant diapirs rise from the base of the lake or from upwelling magma along the lake wall (Fig. 3a). We call this regime the dripping diapirs or simply dripping regime. After $\sim 50$ min, the flow switches to a new regime in which a near-continuous buoyant plume of magma ascends directly from the conduit mouth and rises through the center of the lake (Fig. 3b). The plume oscillates from side to side and accommodates occasional pulses of gas-rich magma. We call this regime the pulsing plume or simply pulsing regime. The time required to reach the latter flow regime motivates the model run time of two hours applied to all parameter variations.

We suggest that the transition between dripping and pulsing regimes is the result of increased mean vesicularity in the lake. Gas-rich magma that reaches the surface is only partially outgassed before being recirculated, resulting in a gradual increase in the lake’s gas content (Fig. 3c). The transition between the regimes occurs when the mean vesicularity is $\sim 0.07$, shortly before the vesicularity reaches a dynamic steady state of $0.085\pm 0.005$ ($2\sigma$). As the mean vesicularity increases, the density contrast between gas-rich upwelling and partially outgassed downwelling magmas decreases. The loss of relative buoyancy reduces the magma ascent rate compared to the conduit inflow and thus suppresses the dripping instability initially present at the conduit mouth. Visual inspection of the maximum and mean velocity time series (Fig. 3d) further suggests that the different surface velocity patterns associated with the two regimes may be a useful metric to discriminate between the regimes in field data.

4.2 Flow regimes

The two flow regimes are controlled by the dimensionless numbers directly affecting the flux of gas into and out of the lake. These are the inflow number, $\mathrm{R_{in}}$ (19), and the outgassing number, $\mathrm{Da_{\phi}}$ (20a), which among other parameters are governed by the inflow vesicularity and rate, and the time scale of outgassing, respectively. We further investigate the physical mechanism underpinning the two regimes and their stability fields by testing the parameter variations given in Table 2 whose results are visualized in the Supplemental material.

4.2.1 Dripping diapirs regime

The dripping diapirs regime is characterized by an intermittent rise of buoyant, gas-rich magma diapirs formed by a dripping instability at the lake inlet. The diapirs either rise directly from the conduit or detach from magma that has ascended as far as one-half to two-thirds of the height of the lake along one of the walls (Fig. 4a-c). Diapirs have a radius of $\sim 1/2$ the conduit radius and rise every $\sim 1.5$–4 min at speeds of $\sim$ 0.2-0.4 m/s.

We find the dripping regime is stable when $\mathrm{Da_{\phi}}$ is high and $\mathrm{R_{in}}$ is low (Fig. 5). $\mathrm{Da_{\phi}}$ is primarily controlled by the imposed surface outgassing time, $\tau_{\phi}$. When $\tau_{\phi}$ is short relative to time scale of convection, gas is removed efficiently from the system and a high density contrast between upwelling and outgassed magma is preserved. While this parameter is particular to our formulation and combines the effects of a range of processes, it highlights the importance of rapid outgassing in preventing gas build-up.

$\mathrm{R_{in}}$ depends linearly on the imposed inflow speed, which remains poorly constrained by field observations. We consider values between 0.1-0.5 m/s based on estimates from Calkins et al. (2008) and Oppenheimer et al. (2009). Alternatively, an increase in inflow vesicularity will slightly reduce $\mathrm{R_{in}}$ since more buoyant magma will rise at a higher ascent speed relative to the inflow. Both slower inflow from the conduit as well as faster convective ascent speed promote stability of the dripping regime.

4.2.2 Pulsing plume regime

In the pulsing regime, buoyant magma rises through a mostly continuous plume of magma near the center of the lake (Fig. 4g-h). The plume does not steadily approach the surface, but oscillates from side to side with time. The plume is occasionally disrupted by diapirs that assemble at the conduit mouth or from the lateral extreme of a bend in the plume. At times, the downwelling flow cuts across the plume near the inlet. Even when continuous, the plume does not maintain a constant flux but accommodates pulses of increased flux ascending along the existing pathway. Pulses arise from a dripping instability similar to diapirs, but with a smaller volumetric imbalance.

We find the pulsing regime is stable when $\mathrm{Da_{\phi}}$ is small and $\mathrm{R_{in}}$ is large (Fig. 5). For slow outgassing, the steady-state lake vesicularity is high, which reduces the ascent speed of gas-rich magma. The regime transition depends on gas in- and out-flux such that the pulsing regime is favored when $\mathrm{Da_{\phi}}$/$\mathrm{R_{in}}\lesssim 0.1$ (Fig. 5). Within the range of outgassing rates we consider, conditions for the pulsing regime are met for $u_{in}\gtrsim$ 0.3 m/s and for $\phi_{0}\gtrsim$ 0.2 We do not test scenarios of $\mathrm{R_{in}}\gtrsim 1$, where inflow exceeds free convection, which may be applicable for advection-dominated lava lakes like Ambrym (Carniel et al., 2003; Lev et al., 2019).

4.2.3 Transitional regime

We observe a transitional regime between dripping and pulsing in which magma ascends as an unsteady plume along the lake walls in the bottom half of the lake before becoming detached and rising through the lake interior to the surface (Fig. 4g-i). The plume is more prone to disruption and diapirs from the conduit are more common in this than in the pulsing regime. Furthermore, the regime is transitional in that it sometimes marks a transition period between dripping and pulsing regimes as gas accumulates in the lake. In the reference case, this transition occurs quickly, in less than a single overturn cycle; however, in other cases it can be more protracted. The transitional behavior can persist until the end of the simulation time for conditions at intermediate $\mathrm{R_{in}}$ (5.0$\times 10^{-2}$ $\lesssim\mathrm{R_{in}}\lesssim$ 7.4$\times 10^{-2}$) and low $\mathrm{Da_{\phi}}$ (Fig. 5). In the cases where transitional behavior is observed after two hours model time, we find the average lake vesicularity continues to increase, suggesting that it may still represent an intermediate stage before the pulsing regime is reached at steady-state.

4.2.4 Stagnant lid

Since the focus of this paper is to understand the surface signatures of different convective regimes in an active lava lake, we are less interested in simulations leading to the formation of a stagnant lid. We observe the formation of a stagnant lid either if surface cooling is high or if outgassing is very inefficient. We find that to avoid a frozen stagnant lid, the cooling time must exceed $\sim$8 hr ($\mathrm{Da_{T}}\lesssim$ 1.1$\times 10^{-5}$). Inefficient outgassing can lead to accumulation of gas-rich magma to the point that the vesicularity of the lake approaches that of the upwelling magma and gas-driven convection ceases (Supplementary Fig. S5e). We consider the limit where bubble segregation through the cooling skin is the only mechanism of gas removal ($\tau_{\phi}=\infty$) and find that process alone is too slow to maintain open convection. In that limit, other processes such as increased permeability through interconnected vesicularity or fissuring of the skin as it deforms must increase the rate of gas removal (Blower, 2001).

5 Comparison to field data

5.1 Surface velocity

We compare simulated surface data to field observations to evaluate which model regimes are applicable to the observed dynamics at Ray Lake. Fig. 6a,c $\&$ e) shows that the dripping and transitional regimes are characterized by more rapid variations in surface velocity than the pulsing regime. Surface velocity field data are calculated using optical velocimetry on thermal imaging observations from December 2012 (see Peters et al. (2014b) and Lev et al. (2019) for details on data collection and velocimetry). Here, we use two representative 1-hr long time series from 14:00-14:59 on Dec 06 and 05:00-05:59 on Dec 26 (Fig. 6i $\&$ j). These sequences were chosen for their low noise, lack of Strombolian eruptions, and to explore a range of behaviors exhibited by the lake.

We identify dominant periods of surface motion using the Fast Fourier Transform of the velocity in both the models and the field data. We exclude lake edges because they do not necessarily participate in the main convective pattern. For model output, we analyze the central 80$\%$ of the lake surface, and for observational data, use an elliptical mask centered on the lake (Fig. 6g). We normalize profiles at each node by their maximum value to reduce the bias toward the center of the lake that experiences higher mean velocities. The absolute values of model velocities range between 0.01 and 0.2 m/s, consistent with previous estimates (Oppenheimer et al., 2009).

The frequency analyses reveal several dominant periods in the simulations that are broadly consistent with those identified in the field data. Simulations in the dripping regime typically show dominant periods around 1.5-4 min, 7-9 min, and 11-15 min (Fig. 7a). The shortest period corresponds to the assembly and ascent of diapirs from the inlet, while the longer periods are related to lateral migration of the upwelling location and lake overturn. Simulations in the pulsing regime lack the shortest period found in the dripping regime but show peaks at 5-9 min, 12-15 min, and occasionally a weaker signal at 18-20 min associated with the lateral migration and disruption of the plume (see Fig. 7b, and Supplementary Fig. S9 for examples with 18-20 periods). Simulations in the transitional regime show periodicities similar to both the dripping and pulsing regimes. Our models predict a periodicity similar to that observed in the field data (Fig. 7e $\&$ f). However, we find that the analysis does not allow a clear distinction of the likely convective regime in Ray Lake. Short-period peaks (1.5–4 min) diagnostic of the dripping regime are masked by noise in the field data and can therefore not be identified reliably enough to settle on a firm interpretation.

5.2 Spreading location

In addition to the velocity magnitudes, we also consider the spatial distribution of the divergence of surface velocity ($\boldsymbol{\nabla}\cdot v$). At each model time step, we track the location of maximum divergence (spreading) at the surface, which we interpret as the upwelling location.

The maximum spreading location time series are shown in Fig. 6b, d $\&$ f. The difference between the regimes is visually apparent: the dripping and transitional regimes are characterized by rapid variations, whereas the pulsing regime shows more continuous migration. We perform an equivalent two-dimensional analysis on the field data by tracking the maximum spreading location over time, in this case omitting the mask used in the velocity analysis to capture the full range of variability (Fig. 6 h $\&$ l). The two observational records show distinctly different behavior. Data from Dec 6 shows greater variability in surface flow directions resulting in spreading locations covering the entire lake surface. The record from Dec 26 shows spreading locations concentrated at a few locations in the lake center or the lake edges arising a single spreading axis in the surface flow field.

Similar to our analysis of surface velocities, we perform a Fourier analysis on time series of the maximum spreading location. The period spectra for field observations are quantified by averaging spectra of the distance from the lake center and of the azimuth relative to an arbitrary reference. Rather than seeking to identify specific short periods diagnostic of dripping, we consider the relative cumulative energy in longer periods ($T\geq 5$ min) of maximum spreading migration as shown in Fig. 8. To reduce the effect of noise, we cut off energy below periods of 1 min, where we do not expect to find variability associated with the convective regime.

As expected based on the velocity magnitude analysis, we find that simulations in the pulsing regime have more energy above $T=5$ min than those in the dripping regime. Within each regime, simulations with lower $\mathrm{R_{in}}$ (slower inflow, faster convection) have more long-period energy. The analysis of observational data returns values in the same range as the simulations, with the two observational records overlapping our predicted regime boundary. We thus interpret that dynamics at Ray Lake appear to straddle the regimes, and that the segment from Dec 6 shows behavior consistent with the dripping regime, whereas the segment from Dec 26 indicates the pulsing regime. Based upon analyses from other dates and times (not shown here), the Dec 26 data seem more representative of the typical behavior of the lake.

5.3 Periodicities from the conduit

Our analysis shows that steady inflow into a convecting lava lake can produce periodic behavior strikingly similar to observations. The question remains whether periodicities in conduit flow remain identifiable in surface observations. To identify how the lava lake modulates imposed conduit periodicities, we test inflow conditions with a sinusoidal time variation in inflow speed of a set period. We test inflow periodicity of 5, 10, and 20 min to cover the observed range and keep all other parameters the same as our reference case. Time-varying conduit inflow result in apparent switching between flow regimes over simulation time. When inflow is slow, dripping behavior is observed, and when inflow is fast, pulsing behavior.

To determine whether the periodic inflow signal is recoverable at the surface, we analyze the surface velocity data in the same way as for constant inflow conditions. Figure 7d shows the surface velocity spectrum recorded for a 5 min conduit inflow periodicity. The imposed periodicity from the conduit is indeed recoverable from the simulated surface velocity. However, the spectrum also shows the excitement of fundamental periods found in the reference simulation at constant inflow. For an imposed 5 min periodicity, we find spectral peaks at 3, 5, and 15 min, while the 7 min periodicity seen in the reference case is suppressed. We conclude that conduit periodicities can be expressed at the surface, but are overlaid with the fundamental periods pertaining to internal lake convection. A clear discrimination of conduit-related from lake-related signals remains challenging.

5.4 Model limitations

Our choice of parameters is informed by field data, but we simplify a variety of processes including the removal of heat and gas from the lake, micro-physical phase interactions, boundary effects, and the potential role of gas slugs. In addition to these processes that are not currently represented in the model, the processes that are included might scale differently in the 2D model than in the 3D lake. Despite these simplifications, we find surprisingly good agreement between our model and field observations across a range of parameters. This may indicate that processes such as diapir formation and flow reorganization may be fundamental to lava lake circulation, and that the predicted behaviour is robust across a range of scales.

With our reference value of $\tau_{T}$, simulations do not reach thermal equilibrium within the set run time of 2 hr. Rather, the mean lake temperature cools at approximately 0.65$\celsius$/hr. If this rate of heat loss continued linearly, we would expect the lake to fully crystallize after $>$120 hr. However, the majority of the heat lost is either by diffusion through the walls where thermal equilibrium cannot be achieved during the simulation time, or at the edges of the lake surface where the layer of magma is thin and rarely disrupted by convection. These top corners are the main source of thermal imbalance, and their contribution will decay with time as they reach thermal equilibrium (Supplementary Fig. S3). Accordingly, the convecting portion of the lake interior does reach an approximately steady-state temperature over the simulation time.

6 Summary $\&$ conclusions

Based on current observational and experimental constraints, our simulations suggest that the dynamics at Ray Lake may straddle the boundary between two convective regimes. In our model, shifts between the dripping and pulsing regimes are triggered by small changes in the relative speeds of conduit inflow to free convection ($\mathrm{R_{in}}$), and rates of surface outgassing to magma transport ($\mathrm{Da_{\phi}}$). Our results suggest that unsteady lava lake convection fed by a constant supply of gas-rich magma from the conduit is capable of generating periodic surface behavior similar to that observed at Ray Lake. We find that surface velocity spectra can be diagnostic of the convective regime, but that the most diagnostic short periods are typically masked by noise in the data. However, spectral energy in the migration of the maximum spreading location at periods longer than 5 min allows to distinguish more reliably between the dripping and pulsing regimes. Our findings further indicate that, using surface velocity data alone, periodic inflow from the conduit cannot be readily discriminated from unsteady lake convection. It is conceivable that observed variations in gas compositions may aid the interpretation of conduit conditions. We conclude that caution should be taken when interpreting surface records at lava lakes as a direct signature of conduit flow.

7 Acknowledgements

This research was supported by the Stanford Earth Summer Undergraduate Program and NSF under award EAR 1348022. T.K. acknowledges support from the Postdoc Mobility Fellowship 177816 by the Swiss National Science Foundation. E.L. was supported by NSF award 1348022.

8 Author contributions

This work is the result of the undergraduate thesis research of J.B., who performed the simulations and data analysis on model output and field data and took the lead in writing and figure preparation with the support of the other authors. T.K. and J.S. developed the research question, T.K. developed the numerical model, and E.L. provided the observational data for Mount Erebus, including image analysis and velocimetry. All authors reviewed and approved the text.

9 Competing interests

The authors have no competing interests to declare.

Supplementary material

We include a table of mineral densities used to compute the average crystal density (Table S1). Functional fit to laboratory measurements of equilibrium crystallinity (Fig. S1). Plot of viscosity as a function of temperature and crystallinity (Fig. S2). Temperature change over a model run showing heat loss through conduction through the lake walls and cooling of the lake corners with constant lake temperature(Fig. S3). Additional parameter variations that modify cooling (Fig. S4), outgassing (Fig S5), and inflow (Fig. S6) conditions and magma rheology (Fig. S7). Resolution testing showing consistent behavior at higher resolution (Fig. S8). Additional selected Fourier analyses in the dripping, transitional, and pulsing regimes and for periodic inflow conditions (Fig. S9).