Numerical simulations of convection in crystal-bearing magmas: A case study of the magmatic system at Erebus, Antarctica

1. Introduction

[2] Characterizing the dynamics of convection currents in magmatic plumbing systems is a challenging task, which has long been approached through experimental [e.g., Jaupart and Brandeis, 1986; Weinstein et al., 1988; Bazarov et al., 2007; Huppert and Hallworth, 2007] and theoretical [e.g., Sparks et al., 1984; Jaupart and Brandeis, 1986; Jellinek and Kerr, 2001] studies that have considered magma as a perfectly homogeneous mixture of melt, crystals and bubbles. More recently, studies based on multiphase theory have been addressing the differential motion of these three phases, taking them in pairs [e.g., Dartevelle, 2004; Dartevelle et al., 2004; Dufek and Bergantz, 2005; Dartevelle and Valentine, 2007; Ruprecht et al., 2008; Dufek and Bachmann, 2010]. These studies have, among other insights, highlighted the complex relationships between convective movements, which tend to keep crystals in suspension because of high melt viscosity, and the gravitational settling of crystals, which drives melt/crystal separation. Most of these numerical and experimental studies, however, have employed highly idealized magmatic systems; only a few have focused on specific natural systems.

[3] Motivated by the fact that they provide a window into the underlying magmatic system, many field studies have been devoted to lava lakes. The nature of convection can be best understood in long-lived systems where low viscosity magmas (e.g., phonolitic, nephelinitic, basaltic andesitic, basaltic) are permanently exposed, such as Erebus volcano in Antarctica [Oppenheimer and Kyle, 2008, and references therein], Nyiragongo in the Democratic Republic of Congo [Demant et al., 1994], Villarrica in Chile [Witter et al., 2004], Erta'Ale in Ethiopia [Bizouard et al., 1980], Pu' u O' in Hawaii [Garcia et al., 2000], and Stromboli in Italy [e.g., Giberti et al., 1992; Laiolo and Cigolini, 2006]. Several aspects of Erebus volcano motivated us to choose it as a subject for a numerical study of convective processes. A unique feature of its lava lake is the large size (up to ∼10 cm) of crystals it contains [Kelly et al., 2008]. These megacrysts represent more than 97 vol.% of the crystal cargo and beg the question as to whether their sedimentation interferes with convection or if they act as passive tracers with trajectories diagnostic of a particular convective regime. Many studies have been devoted to Erebus, providing a wealth of data on which numerical simulations can be based. They provide indications of lake and upper conduit dimensions [e.g., Dunbar et al., 1994; Aster et al., 2008; Dibble et al., 2008; Kyle et al., 1992], the current low rate of crystallization [Dunbar et al., 1994; Kelly et al., 2008], as well as good constraints on the magmatic properties and characteristics such as crystal content and size [e.g., Kelly et al., 2008; Dunbar et al., 1994], temperature [e.g., Calkins et al., 2008; Kelly et al., 2008; Sweeney et al., 2008], viscosity [Sweeney et al., 2008], crystal age [e.g., Reagan et al., 1992] and steadiness of convection [e.g., Kelly et al., 2008].

[4] Here, we reduce the gap in numerical studies dedicated to a given volcano by studying the behavior of the shallow magmatic system of Erebus (from reservoir to lava lake) when it is assumed to contain either a single-phase mixture of silicate melt and crystals, or a bi-phase combination of silicate melt plus crystals. Carrying out such numerical simulations involves two complementary aspects. From a modeling perspective, contrasting the convective motions generated by a simplified mixture to those generated by a bi-phase system (where crystals can settle) identifies which features of magmatic convection are best viewed using a single- or a multiphase system. This represents an important step in the development of numerical models that include the three phases (melt, crystals, and bubbles) present in magmas. The simulations also allow us to tackle technical issues concerning boundary conditions, which govern the long-term behavior of the system. More specifically, our study is relevant to understanding the behavior and stability of Erebus lava lake. A sustained activity is surprising for a lake of such small dimensions; it implies the existence of a vigorous convecting system that maintains the lake at a quasi-constant temperature. Considering the fluidal nature of molten phonolite, it is reasonable to ask whether temperature gradients alone can generate sufficient buoyancy to keep the lake from freezing, or if additional buoyancy sources, such as gas bubbles, need to be involved. The underlying driving mechanism is not unique to Erebus, and its understanding might shed light on other volcanoes.

[5] In the next section, we briefly describe the bi-phase model used and basic assumptions underlying its application to magmas. 6 concerns validation and verifications tests of the single-phase part of the model via reproduction of analogue experiments of transient convection, and of the bi-phase part by reproduction of the bulk behavior of a simple Poiseuille flow. In 11 we describe the different numerical configurations tested to study the Erebus magmatic system and the results of the numerical simulations. Overall, these results suggest the existence of a minimum conduit diameter that ensures effective and lasting convection, and reveal partial settling of crystals - provided they are considered as a separate phase - leading to the formation of an enriched layer of crystals at the bottom of the magma chamber. In 25 we compare our results with other theoretical and experimental works to draw conclusions on the suitability of using multiphase physics to simulate magma convection. We find that the extension of our results to Erebus and other magmatic systems requires consideration of bubbles.

2. Physical Model Description

[6] The dynamics of crystal-bearing magmas can be modeled as bi-phase flows (a mixture of crystals/solids/dispersed phase and fluid/continuous/carrier phase) using either a continuum (Eulerian) or a discrete-particle (Lagrangian) approach. Here, we will consider a system consisting of phases sharing the same pressure field (Eulerian-Eulerian approach) [Dartevelle, 2004; Dartevelle et al., 2004]. We wished to limit the complexity of the model so as to capture the fundamental differences between mixture theory, where crystals and melt are considered as a single phase with modified properties, and the bi-phase formulation that considers crystals as a separate phase. This led us to make two fundamental assumptions, as follows.

[7] First, we considered Erebus magma as a bi-phase flow containing melt and crystals, thereby ignoring the presence of a gas phase. This is a shortcoming because the sustained gas emission from the lake surface and sporadic Strombolian explosions through the lake are clear indications that gas bubbles play an important role at Erebus [e.g., Giggenbach et al., 1973]. Ignoring gas will lead to underestimation of the convective dynamics of the simulated system, slowing down motions because of the neglected buoyancy source. While this severely limits the applicability of our model to the real system in terms of making direct comparisons between model outputs and observed parameters, we nevertheless can simulate a magmatic system with many other realistic characteristics (complex geometry, physical parameters of the magma, timescales, etc.) and draw conclusions on the role of megacrysts in magmatic convection.

[8] Second, we treated crystals as rigid spherical particles of constant size in order to realize consistent rheological formulations for the mixed and bi-phase systems. This is a significant shape simplification because Erebus megacrysts have aspect ratios (ratio between long and short axes) of 2–3 [Dunbar et al., 1994], an anisometry sufficient to affect the rheology of the suspension [e.g., Mueller et al., 2010]. The works of Caricchi et al. [2007] and Costa et al. [2009] have treated magma as a homogeneous mixture of crystals and melt. They proposed effective viscosity formulations taking into account the non-sphericity of the particles through the linear dependence between the average aspect ratio and maximum packing fraction given by Pabst et al. [2006]. Matching such dependence with a bi-phase flow, however, is beyond the scope of the present work, mainly because drag laws and solid moment tensors for describing the behavior of non-spherical crystals are not currently available in the literature. It can nevertheless be estimated that considering our crystals as rigid spheres underestimates the viscosity of the suspension [Picard et al., 2010].

[9] The numerical code that we used to solve the bi-phase flow equations is MFIX Ver. 2.0, 2004 (Multiphase Flow with Interphase eXchanges), a program developed at the Department of Energy National Energy Technology Laboratory of U.S. The software package VISIT Ver. 1.12 developed by the U.S Lawrence Livermore National Laboratory was used to visualize simulation outputs and all numerical data were post processed by script shells developed in MATLAB^® and FORTRAN. MFIX is based on several works developed since 1967 [Syamlal et al., 1993; Syamlal, 1994, 1998]. The code, based on a discretization using staggered and finite volume methods, has provided a platform by which to explore various volcanic phenomena: (i) as a closed and isothermal system to simulate magma mixing [Ruprecht et al., 2008], (ii) as an open and isothermal system to explore conduit flow [Dufek and Bergantz, 2005], (iii) as an open and non-isothermal system to explore the magma chamber dynamics [Dufek and Bachmann, 2010] and conduit flow [Dartevelle and Valentine, 2007]. In this study we used MFIX to study the magmatic system feeding the lava lake at Erebus volcano as (i) closed and non-isothermal, and (ii) open non-isothermal systems. This is the first study using a closed and non-isothermal system. Since the temperature gradients in the simulated magmatic system are very large, this approach has limitations because the overall cooling causes melt contraction, which is incompatible with the finite-volume numerical scheme applied to a closed system. This issue is addressed in 20.

2.1. Governing Equations

[10] MFIX uses the Navier-Stokes equations for momentum together with supplemental equations for the interphase forces, plus mass conservation and energy equations, which enable the simulation of incompressible multiphase flow with a Newtonian carrier phase. The main transport vector equations for a bi-phase system without chemical reaction between phases are reported in Table 1. The constitutive terms of each equation are explained in the text and the nomenclature is summarized in Table A1. A positive sign convention is adopted throughout this study, which means that stresses are positive when both the outward normal and force act in the same sense relative to the coordinate system x and y, or East and up, respectively.

Table 1. General Transport Equations of the Physical Model

Equation	Number
Continuity
εm + εs = 1	(14)
Melt
	(15)
Solid
	(16)
Momentum
Melt
	(17)
Solid
	(18)
Stress tensor
Melt
	(19)
Granular
+ τsP if εm ≤ εm* Plastic/Frictional	(20)
= −Psv  + τsv if εm > εm* Kinetic - Collisional	(21)
Granular normal stress
PsP = 1024(εm* − εm)10εs	(22)
Psv = 2(1 + e)ρsgoεsΘs	(23)
where go = f(εm, εs)	(24)
Θs = f(K1, K2, K3, K4, εs, tr( ))	(25)
K1 = f(e, go, ρs)	(26)
K2 = f(e, go, ρs, dp, εs, K3)	(27)
K3 = f(e, go, ρs, dp, εs)	(28)
K4 = f(e, go, ρs, dp)	(29)
Melt viscous stress
	(30)
Granular viscous stress
	(31)
	(32)
	(33)
Momentum interface transfer coefficient
Drag forces
Fms = f(εs, εm, ρm, dp, Res, Vr, vs, vm)	(34)
Terminal velocity
Vr = f(εm, Res)	(35)
Particle Reynolds Number
Res = f(dp, ρm, μm, vs, vm)	(36)
Energy
Melt
	(37)
Solid
	(38)
Heat conductivity
Melt conductivity
qm = εmkm∇Tm	(39)
Granular conductivity
qs = εsks∇Ts	(40)
Heat interface transfer coefficient
γms = f(km, εs, Nu, dp)	(41)
Nusselt Number
Nu = f(εs, εm, Res, Pr)	(42)
Prandtl Number
Pr = f(Cpm, μm, km)	(43)

[11] Continuity requires that the sum of the volume fractions always equals one (equation (14); equations (14)–(43) are given in Table 1). The differential form of the mass conservation equation (equations (15) and (16)) states that the accumulation of mass per unit volume plus the convective flux of mass per unit volume (first and second terms of LHS) is equal to zero (RHS). This absence of phase exchange means that there is no crystal growth or nucleation. The Navier-Stokes equations account for the rate of increase of momentum per unit volume and for the rate of momentum transferred by convection per unit volume (first and second terms of LHS in equations (17) and (18)). Thus the whole LHS term is equal to (in the RHS from left to right) (i) the interaction within the phase accounted for as the stress tensor (normal and shear stress components) per unit volume, (ii) the net gravitational body force acting on each phase, and (iii) and (iv) the interaction forces between the fluid and solid phases per unit volume (drag and buoyancy force). The term (iv) only exists in equation (18). The formulation of the stress tensor for the melt (equation (19)) follows a Newtonian form. The granular phase (equations (20) and (21)) follows the transitional theory from viscous to plastic regime proposed by Johnson and Jackson [1987] through the void fraction at minimum fluidization ε_m*. Thus, at high concentrations, the grains impinge on each other, losing mobility and giving rise to a frictional (plastic) dissipation and transfer of momentum (equation (20)) [Schaeffer, 1987]. Otherwise (equation (21)), the grains either move randomly by kinetic transport (particle translation) or by grain-grain collisions following the Chapman-Enskog approach for dense gases [Chapman and Cowling, 1970]. Grains interact with the fluid through drag laws [Jenkins and Savage, 1983; Lun et al., 1984; Ding and Gidaspow, 1990]. Equations (22)–(36) define the factors composing the stress tensor for the melt and granular phase, and the drag forces [Syamlal et al., 1993]. The LHS of energy equations (37) and (38) indicates that the net rate of change of temperature within a grid cell is equal to the sum of the work resulting from the RHS: (i) heat conduction within the phase and (ii) the heat transfer between phases. In this study we have neglected the viscous dissipation and the work due to drag forces. Equations (39)–(43) define the factors relevant to (i) and (ii) [Syamlal et al., 1993].

[12] In simulations using two phases, crystals have a constant density and the melt has a temperature-controlled density. In single-phase simulations, the mixture phase density, ρ_(m+s), is modeled as the weighted sum of the constant density of the dispersed phase and the thermally controlled density of the continuous phase:

Where ε_m is the liquid fraction, ρ_s is the crystal density, T_m is the temperature of the melt, T_o is the initial temperature, ρ_o and α are the density and thermal expansion coefficient of the melt at initial temperature.

[13] The crystal content may play a fundamental role in increasing the magma bulk viscosity, μ_(m+s). Close to the maximum packing fraction (1 − ε_m*), the effect of crystals dominates over the liquid phase, the suspension loses mobility and the bulk mixture viscosity increases significantly. Among others, Caricchi et al. [2007] determined that the presence of crystals causes the mixture to react as a non-Newtonian fluid. This behavior could result in visco-plastic or Bingham-plastic behavior, where the magma exhibits a finite yield strength [e.g., Murase and McBirney, 1973]. However, the modeling of this behavior remains uncertain and we have instead adopted the broadly accepted rheology for the mixture given by Krieger and Dougherty [1959]:

Where μ_m is the dynamic viscosity of the melt and ε_s is the crystal fraction. Since we consider mono-dispersed (same size) crystals as rigid spheres without attractive or repulsive forces, the Einstein coefficient was fixed at the theoretical value, [η] = 2.5. We note that the more sophisticated effective viscosity formulations proposed by Caricchi et al. [2007] and Costa et al. [2009] converge on the Krieger and Dougherty [1959] relationship when crystal fraction is lower than the maximum packing fraction. As mentioned in 6 and 11, our simulations scarcely reach such crystal fraction.

2.2. Thermal Boundary Conditions

[14] To define the heat that flows out of the computation domain, MFIX has two kinds of wall boundary conditions for the energy equations: (i) non-conducting wall (no interaction with the surrounding medium) in which the wall temperature, T_wall, is fixed, and (ii) boundary heat flux specified by a heat transfer coefficient C_m:

Where ∂/∂_n denotes differentiation along the outward-drawn normal, and is the wall temperature as a function of distance.

[15] Algebraic methods can be used for the calculation of C_m across plane walls. To calculate how much heat is lost through the walls, we have used the algebraic model of Turcotte and Schubert [2002], which states that the temperature of the wall rock as a function of time can be obtained by solving the one-dimensional, time-dependent heat conduction equation for a semi-infinite half-space, initially at uniform temperature and suddenly brought to a new temperature that remains constant from then on. The rate of heat transfer depends on the temperature gradient and the thermal conductivity of the material as follows:

Where T_o is the initial temperature of the system (the same value as that used in equation (1)), T_fwall is the temperature far away from the wall along the outward-drawn normal, x is the distance from the wall, t is time and κ_wall is the thermal diffusivity of the wall rock. At long times, the system becomes insulated because the rate of heat transfer slows down. We calculated the derivative of equation (4) to obtain C_m.

2.3. Numerical Considerations

[16] Optimal grid size can be estimated using scaling analysis. When the local Rayleigh number reaches a critical value (Ra_c) [e.g., Sparrow et al., 1970; Jaupart et al., 1984], a convective instability develops and the thermal boundary layer breaks down into incipient plumes. The theoretical length of the grid size of the cell that captures well the dimensions of such a plume is given by considering the Rayleigh number, Ra, at its critical value (Ra = Ra_c):

Where ρ_m is the melt density, ΔT is the driving temperature difference, g is the gravity, κ_m is the thermal diffusivity of the melt, μ_m is the dynamic viscosity and dy is a characteristic length (in this case cell size).

[17] Several computing parameters control numerical convergence. In MFIX, time progression is governed by three parameters: starting, maximum and minimum time steps. MFIX uses an automatic step adjustment, which can be upwards to the maximum time step or downward to the minimum time step, which ensures that convergence is obtained [Syamlal, 1998]. The minimum time step has a physical limit. At the sub-grid scale, in an individual computational cell, particle acceleration within a given time step should be uniform. Sub-grid crystal velocity is described by the Newtonian equation of motion:

Where V_l is the instantaneous particle velocity, v_s is the Stokes velocity, M is the mass of the particle, ξ is a proportionality constant related to the drag forces exerted by the melt (ξ = 3πd_pμ_m), which depends on the particle diameter d_p. The minimum time step so that the solid phase speed is a macroscopic property is defined as the time taken for steady motion to be reached, which is when crystal velocity equals 99.5% of its steady Stokes settling velocity, defined as v_s = (Δρgd_p²)/(18μ_m), where Δρ is the density difference between particle and fluid.

3. Validation and Verification of the Physical Model

[18] When modeling convective phenomena, dynamical parameters such as cooling rate, the wavelength of instabilities that form during overturn, the velocity of the resulting plumes, and the evolution of the bulk viscosity of the mixture are key criteria to verify the numerical code. For single-phase simulations, we validated our model against an analogue experiment performed by Jaupart and Brandeis [1986]; this validation yielded 2D numerical solutions of fluid convection. We used this model outputs to verify that descending plumes have velocities in conformity with the formulation proposed by Read [1984] and Youngs [1984]. The wavelength of instabilities can be predicted, as suggested by Odé [1966] and Lister and Kerr [1989]. To verify our code for bi-phase simulations, we used a fully developed Poiseuille flow to retrieve a bulk viscosity that we compared to the mixture formulation proposed by Krieger and Dougherty [1959] in equation (2). We detail validation and verifications below.

3.1. Validation and Verification of Momentum and Heat Transfer for Single-Phase Simulations

[19] For single-phase simulations, we tested whether MFIX could reproduce data, observations and conclusions of previously reported analogue experiments of Jaupart et al. [1984] and Jaupart and Brandeis [1986] that were aimed at simulating convection in magma chambers. The analogue experiments reproduced thermal convection by imposing a thermal gradient in a tank full of various silicone oils. We focused on the experimental convection of three kinds of oil, the physical properties of which are reported in Table 2a.

Table 2a. Parameters of the Experiments^a

Parameter	Symbol	Unit	47V20	47V20b	47V500	47V1000
Initial temperature	To	°C	49.4	49.4	54.4	36.5
Temperature of top and bottom walls (BC1 and BC4)	Tf	°C	27.6	27.6	27.2	21.3
Temperature contrast (To − Tf )	ΔT	°C	21.8	21.8	22.2	21.8
Density	ρm	kg/m3	930	930	942	946
Dynamic viscosity	μm	Pa s	0.012	0.12	0.284	0.767
Heat capacity	Cpm	kcal/ kg °C	0.368	0.368	0.371	0.368
Thermal conductivity	km	W/m °C	0.145	0.145	0.144	0.145
Thermal diffusivity	κm	m2/ s	1.01 × 10−7	1.01 × 10−7	9.87 × 10−8	9.96 × 10−8
Thermal expansion	α	1/°C	1.07 × 10−3	1.07 × 10−3	9.45 × 10−4	9.45 × 10−4
Rayleigh number	Rac	dimensionless	1.7 × 108	1.7 × 107	6.9 × 106	2.5 × 106
Observed height of the upper thermal boundary layer	hc	m	0.21 × 10−2	0.48 × 10−2	0.60 × 10−2	0.87 × 10−2
Theoretical height of the upper thermal boundary layer	dyd	m	0.21 × 10−2	0.45 × 10−2	0.61 × 10−2	0.86 × 10−2
Error % for the height of the upper thermal boundary layer			0	6	2	1
Observed wavelength of the cold crust	λobs	m	0.60 × 10−2	1.50 × 10−2	1.80 × 10−2	2.70 × 10−2
Theoretical wavelength	λce	m	0.58 × 10−2	1.48 × 10−2	1.77 × 10−2	2.57 × 10−2
Error % for wavelength			3	6	2	5
Number of starting plumes as a function of viscosity		dimensionless	32	15	10	6
Theoretical number of starting plumes		dimensionless	34	15	12	8

• a The data are based on Jaupart and Brandeis [1986] and Bluestar Silicones France SAS (http://www.bluestarsilicones.com/).
• b 47V20 oil with modified dynamic viscosity.
• c Value calculated with equation (5) and dy is the thickness of the tank.
• d Value calculated with equation (5) and Ra = Ra_c = 1600 [Jaupart et al., 1984].
• e Value calculated with equation (9) from Odé [1966] and Lister and Kerr [1989].

[20] In order to simulate conditions of instantaneous cooling (i.e., a temperature drop ΔT at time t = 0), the analogue experiments started with a hot and isothermal fluid layer at T_o = 49°C enclosed in a Plexiglas tank capped above and below by copper plates. A cryostat maintained the cooling fluid at 4°C in closed circuit. At time t = 0, both copper plates surrounding the hot fluid tank were linked to the cryostat circuit. After about 3 min, the temperature of both copper plates had dropped by a certain value. By setting the temperature of the cryostat close to that new value, Jaupart and coworkers were able to stabilize the coolant temperature to within 0.1°C in less than 6 min. The temperature of both upper and lower boundaries was thus maintained at a constant value T_f = 27.6°C for the rest of the experiment. Results show that the cooling rate is important over the first 3 min, before the temperature of the cryostat is raised to a value close to the final temperature of the tank. The cooling rate then diminishes significantly.

[21] The fluid numerical simulation domain consists of a 2D rectangular system with x- and y- origins at the bottom left-hand corner of the tank (Figure 1a). The aspect ratio of the experimental tank was preserved, although the simulations were carried out in 2D. All the boundaries were treated as “non-slip walls” (NSW). For the simulated range of temperatures (27.6°C to 49°C), we found that the size of the cell required to resolve the time at which instabilities begin to grow is around 0.01 cm (equation (5)). We searched for the optimal computational grid size by trial and error accordingly, between 0.125 and 0.01 cm. As expected, we observed that the plumes (convective instabilities) were better defined as grid size was decreased. The resulting average temperature in the central part of the tank is constant (within ∼3°C) for the three simulations with grid sizes between 0.01 and 0.03 cm. This is close to the theoretical value, and we selected a grid cell of 0.03 cm (960 × 512 computational cells) for better computational efficiency. A satisfactory match between experiments and runs will be taken as the sign that 2D simulations of a 3D tank with a square upper surface can capture important features of the convective process. These features include the vertical temperature profile, the evolution of the temperature of the well-mixed central part of the tank, the progression of the unstable cold front, and the wavelength of the instabilities that compose it.

3.1.1. Comparison Between Numerical and Experimental Results

[22] We first describe the simulation with a tank full of oil 47V20 (Table 2a and Figure 1). The initial oil temperature was set at a uniform and constant value of 49.4°C. The tank was then cooled from above and below, which means that BC1 and BC4 were fixed at a constant temperature of 27.6°C while the vertical walls BC2 and BC3 were considered as insulators (C_m = −10⁻⁸°C/m, Figure 1a). Following Jaupart and Brandeis [1986], we compare numerical runs and experiments as a function of a dimensionless time t* = t/τ, with τ as a diffusive thermal scale, and a dimensionless temperature T* = 1 − (T_o − T_m)/(T_o − T_f), which is a function of the initial fluid temperature, T_o, the horizontal average fluid temperature at a given depth, T_m, and the plates' temperature, T_f.

[23] As shown in Figure 1b, the average temperature profiles from the analogue experiments have steeper profiles at the boundaries than those of the numerical runs (especially near the top of the tank). The temperatures of the central part of the tank, in the well-mixed layer that is nearly isothermal, are, however, identical for both the model at t* = 0.21, 0.39 and 0.63 and the experiment at t* = 0.29, 0.5 and 0.79. These dimensionless times correspond to approximately 9, 16, and 27 min, and 12, 21, and 33 min for model and experiment, respectively. This yields an initial time shift of 3 min, which corresponds well with the time taken by the system to stabilize before the cryostat is set to the temperature of equilibrium. Jaupart and Brandeis [1986] calculated a theoretical dimensionless function T_m* for the variation in temperature of the well-mixed layer (the constant flat profile for each curve in Figure 1b). This dimensionless function is plotted in Figure 1c and compared to both analogue and numerical simulations. We observe a discrepancy between the three curves. The discrepancy between theory and analogue model was already noticed by Jaupart and Brandeis [1986]. They report that experimental data fall on the theoretical curve when shifted in time because the theoretical curve assumes instantaneous cooling of the boundaries. The same applies to our simulation (“numerical translated” curve in Figure 1c), and explains why our simulation approaches the theoretical curve at earlier times.

[24] Rayleigh-Taylor (RT) instabilities occur in a fluid cooled from above when the surface reaches a temperature threshold that generates sufficient density contrast to induce convection. The time sequence of the development of the convection for the experiment and simulation at early times is shown in Figure 1d. The general shape and sizes of the instabilities match quite well, with a time lag varying from 3 min to 6 min at later times. The more numerous instabilities in the analogue experiments can easily be explained by the fact that, unlike our 2D simulation, the shadowgraphs capture several layers of instabilities along the line of sight. Read [1984] and Youngs [1984] found that the cold front formed by the plumes sinks into the warmer fluid as:

Where h(t) is the height of the plume envelope as a function of time, t, β is a constant equals to 0.05 [Dimonte, 1999; Dimonte and Schneider, 2000], g is the gravity constant (taken as 9.8 m/s²) and A_T is the Atwood number defined as:

Where ρ_cold and ρ_hot are the densities of the plumes and warm fluid surrounding them, respectively. We recorded in our simulations the downward progression of the dense plumes and the depth and time at which they vanish, when the density contrast becomes negligible. We found that plume velocity follows a slightly parabolic profile (R² = 0.99) (Figure A1a), which at late time approximates a straight line whose slope represents an average speed of 0.3 cm/s. We verified that this profile approaches that given by the theoretical growth predicted by equation (7) for the simple case where the Atwood number is considered constant with ρ_cold = 0.95 at T_f = 27.6°C and ρ_hot = 0.93 at T_o = 49.4°C. Jaupart and Brandeis [1986] report a constant plume velocity (Figure A1b) in the order of 0.5 cm/s at early times, which matches our results within experimental uncertainty. They noticed that the plumes will become slower as cooling proceeds to reach 0.3 cm/s in the region of the well-mixed layer and as they approach the lower boundary layer.

[25] We increased the viscosity of the oil 47V20 by a factor ten and ran a simulation identical to the one described above (Table 2a). We also simulated other oils (oil 47V500 and 47V1000), which led us, like Jaupart and Brandeis [1986], to vary the temperatures of the initial and boundary conditions BC1 and BC4. We observe that, for high viscosity values, RT instabilities look like fingers rather than ripples (compare Figure A2a with Figure A2d). The number of instabilities developing at very early times decreases as viscosity increases (Table 2a), following a logarithmic profile (Figure A2e). We conclude that, for given initial and boundary conditions, the number of developing instabilities is roughly inversely proportional to bulk viscosity. This is confirmed by measurements of the wavelength, λ_obs, that this conductive layer forms at the onset of the first instabilities (e.g., Figure A1f). Odé [1966] and Lister and Kerr [1989] used scaling analysis to quantify such wavelengths, λ_c:

Where μ_hot and μ_cold are the viscosities of the hot and cold fluid, and h_c is the thickness of the conductive layer (e.g., Figure A1g). For μ_hot ≪ μ_cold, which is always the case here, this theory defines the value of the function Q(μ_hot/μ_cold) = 0.94. The observed and theoretical values of the wavelengths match within a margin of 4% (Table 2a). Knowing the ratio between the width of our domain and λ_c, we can estimate the number of theoretical instabilities that a simulation should develop. This number is close to the one that we observe (Table 2a).

[26] From our data, we conclude that the theoretical thickness, dy, of the upper boundary layer calculated through equation (5), when Ra = Ra_c =1600 [Jaupart et al., 1984] yields values that match our results within a margin of 2% (Table 2a). This calculation provides anotherverification of our 2D model and confirms the accuracy of our simulations within the complexity of the 3D experiments and theoretical values (e.g., Ra_c).

3.2. Verification of Momentum Coupling Between Melt and Solid for Bi-phase Simulations

[27] We verified bi-phase simulations by calculating the equivalent bulk viscosity of melts carrying various amounts of particles and comparing it to the one obtained through the mixture theory of Krieger and Dougherty [1959]. The equivalent bulk viscosities were obtained by simulating a plane parallel Poiseuille flow with varying crystal loads.

[28] We defined a 2D horizontal domain of 10 × 0.5 cm² containing an arbitrary silicate melt enriched in small particles representing crystals at its center (Figure 2a). The simulations were isothermal at 900°C and the physical properties of crystals and melt are listed in Table 2b. The initial conditions were composed of a 8 × 0.5 cm² crystal-bearing center surrounded by two 1 × 0.5 cm² crystal-free regions that ensured that the solution was not perturbed by boundary conditions (Figure 2a). Null velocity was set for both phases and the entire domain was initially at hydrostatic pressure along the y-direction. BC1 and BC2 were set as NSW, causing the velocities of melt and crystals to approach zero on the horizontal boundaries, and BC3 and BC4 were set as fixed-pressure inlet and outlet, respectively. We defined a grid size of 0.01 × 0.01 cm² (1000 × 50 computational cells).

Table 2b. Parameters Used in the Plane-Parallel Poiseuille Flow With a Particle-Laden Suspension^a

Parameter	Unit	Crystals		Melt
		Symbol	Value	Symbol	Value
Density	kg/m3	ρs	2.508 × 103	ρo	2.19 × 103
Dynamic viscosity	kg/m s	n/a	n/a	μm	1.38 × 104
Volume fraction	dimensionless	εs	0–0.6	εm	0.4–1.0
Packed bed (minimum) void fraction	dimensionless	n/a	n/a	εm*	0.05
Maximum packing fraction	dimensionless	n/a	0.95	n/a	n/a
Particle diameter	m	dp	3.74 × 10−5	n/a	n/a
Restitution coefficient between particles	dimensionless	e	0.8	n/a	n/a
Angle of internal friction	° (degree)	ϕ	0	n/a	n/a

• a Labels n/a mean not applicable.

[29] We performed seven simulations where the only difference was the crystal volume fraction, which varied from 0 to 60 vol.% in steps of 10 vol.%. The simulations develop a Poiseuille flow between parallel plates in response to the pressure gradient between BC3 to BC4. For simulations containing from 0 to 30 vol.% of crystals, the inflow boundary at BC3 was set at constant pressure of 2 × 10⁴ Pa and the outflow boundary at BC4 to 100 Pa. For the simulations containing 40 to 60% particles, decreasing the pressure gradient was necessary to reach convergence; otherwise some cells in which solids were loosely packed reached large pressure and velocity values, resulting in unstable calculations. Thus for 40, 50, and 60 vol.%, we fixed a pressure inlet of 2.5 × 10³, 200, and 100 Pa, and an outlet pressure of 100, 100 and 60 Pa, respectively.

3.2.1. Comparison Between Steady State, Bi-phase Flow and Mixture Theory

[30] In steady state, the bulk behavior of the mixture can be considered because particles and fluid have nearly the same velocity. We tracked the horizontal velocity of a vertically centered cell (at coordinate 8; 0.25) to determine when our simulations reach steady state because the time evolution of velocity of this cell is constant in steady state (see the gray-shaded area in Figure 2b). In the x-direction, the governing Navier-Stokes equation for such a mixture reduces to:

Where P_(m+s) and U_(m+s) are, respectively, the pressure and horizontal velocity of the bi-phase flow at every cell along the y-direction. To calculate the bulk dynamic viscosity, μ_(m+s), we determined the pressure gradient (the first term of equation (10)) between the two consecutive faces of a range of cells at coordinates [8; y] and [8 + dx; y]. Then, we determined the double derivative of the velocity in the y-direction (the second term of equation (10)) along coordinates [8; y]. There was no single time at which the value of the bulk viscosity of the bi-phase flow could be determined for all seven simulations. We thus decided to perform this evaluation after steady conditions have been reached, but a little before the bi-phase flow leaves the bin (e.g., Figure 2c for 20 vol.% crystals). This yielded bulk viscosities at 500, 1500, 1500, 1400, 1500, 2000, 1200 s for 0, 10, 20, 30, 40, 50 and 60 vol.% of crystals, respectively. We repeated the calculations for each of the above times plus 50 and 100 s to estimate the error of our measurements. The resulting bulk viscosities are shown in Figure 2d. The error bars in Figure 2d represent the ranges of minimum and maximum values of our three measurements for a single simulation, while the data point is their average value, which we compare to the viscosity defined by Krieger and Dougherty [1959] (equation (2) with ε_m* = 0.05). In these simulations, crystal segregation leads to significant changes in the pressure gradient and the velocity field. As a result, the bulk viscosities calculated with equation (10) only match approximately those of the perfect mixture assumed by equation (2).

4. Erebus Lava Lake Simulations

4.1. Physical Parameters of the Natural System

[31] Erebus volcano on Ross Island is the southernmost active volcano in the world (Figure 3a). The edifice has a wide base (33 km in diameter) with gentle slopes (13°) that are truncated at an elevation of approximately 3700 m by a 600-m-wide summit crater. Inside lies the 200 m diameter Inner Crater that hosts the world's only active phonolite lava lake, informally called Ray's lake in the recent literature, which has been active since at least 1972. Another small lake at the site of Werner's Vent was briefly present in 2004, while another vent called Active Vent has been the site of occasional small ash-producing explosions [Calkins et al., 2008]. Our study concentrates only on the persistent Ray lava lake (Figure 3b). Examination of the surface elevations of Ray Lava Lake, Werner Vent and Active Vent indicates that they were more or less at the same altitude in December 2001 [Csatho et al., 2008].

[32] An extensive series of studies have been carried out on Erebus [Oppenheimer and Kyle, 2008, and references therein]. Among them, geochemical and mineralogical studies were devoted to the phonolite lava. Kelly et al. [2008] studied lava bombs erupted between 1972 and 2004, and the lavas erupted from the summit region in the last 17 ka. They found that the mean composition of matrix glass from lava bombs is characterized by 55.35 wt.% of SiO₂, 9.04 wt.% of Na₂O, 5.64 wt.% of K₂O and 1.87 wt.% of CaO. Lava bombs and summit lava textures are similar. Phonolite lavas from summit region have ∼30% of anorthoclase feldspar crystals with sizes up to ∼6 cm in length, while lava bombs have 30–40 vol.% of anorthoclase crystals with sizes up to ∼10 cm in length, both along with minor quantities (<3 vol.%) of other minerals (olivine, opaque oxides, clinopyroxene and apatite). The mean compositional range of the anorthoclase crystals is An_16.2Ab_65.8Or_17.9 (An = Anorthite; Ab = Albite and Or = Orthoclase). No systematic study of the aspect ratios of the crystals at Erebus has been carried out; however the average length is commonly twice to 3 times the width [Dunbar et al., 1994]. Melt inclusions trapped in the crystals record homogenization temperatures around 1000°C [Dunbar et al., 1994]. The water contents of these anorthoclase-hosted inclusions range from 0.12 to 0.39 wt.% [Seaman et al., 2006]. The Erebus phonolite density is estimated to 2400 kg/m³ [Eschenbacher, 1998].

[33] Observations combined with geophysical, geochemical, and petrological studies provide some constraints on the geometry and dimensions of the magma system and convective process itself. For example, Csatho et al. [2008] used airborne laser scanning (lidar) and Calkins et al. [2008] used thermal infrared images to estimate the lake surface area as ∼770–1400 m² and its diameter as ∼30–40 m (Figure 3) in the period from around 2001–2004. However, the lake level does vary over time with corresponding changes in surface area [Calkins et al., 2008]. In 2009, for instance, the lake area increased to 2500 m². Field observations and videos of the lava lake revealed that hot lava reaching the surface of the lake was rapidly cooled on contact with air, resulting in a thin crust beneath which molten lava kept convecting [Oppenheimer et al., 2009]. Kyle et al. [1992] suggested that the continuous activity at this large alkaline phonolitic volcano implies a deep source of melt rising through a narrow vertical conduit. Calkins et al. [2008] based on the average temperature of the lava lake estimated that a minimum magma flux of 150 kg/s is needed to sustain the thermal output of the lake. Then by assuming a simple Poiseuille flow, they calculated that ascending magma requires a minimum conduit diameter of ∼4 m. Dibble et al. [2008] and Oppenheimer et al. [2009] reported a funnel shaped crater ∼20-m-deep after Strombolian explosion evacuated the lake, tapering to a single vent (uppermost conduit) of 5–10 m diameter. Aster et al. [2008] found that very long period seismic signals accompanying Strombolian explosions are related to the gas slug ascent in the conduit and are consistent with a kink in the magmatic system at around 400 m beneath and horizontally offset from the lava lake. A study focused on the melt inclusions trapped into the crystals revealed that they grew between depths as shallow as 400 m and the surface [Dunbar et al., 1994], though more recent modeling of the melt inclusion CO₂ and water contents suggested entrapment depths corresponding to pressures of up to 300 MPa [Oppenheimer et al., 2011]. Sweeney et al. [2008] found that, for a crystal content of 30 vol.%, the radiated heat flux associated to the SO₂ emission rate corresponds to a cooling of 65°C between upwelling and downwelling magma and is sufficient to drive convection. This temperature difference between magma entering and leaving the lake is consistent with that reported by Calkins et al. [2008] who estimate that ΔT must be inferior to 120°C to sustain convection for same crystal content. Sweeney et al. [2008] also estimated that, as pressure decreases, viscosity varies significantly (from 10⁴–10⁹ Pa s) as a function of water and above all crystal content.

[34] In addition to the studies mentioned above, mineralogical and isotopic composition of lava and tephra suggest that the Erebus magmatic system reached a steady state or is close to its equilibrium, with low to insignificant crystallization taking place. For example, Calkins et al. [2008] report a steady lake level over periods of weeks, constant degassing and stable mean surface temperature. Kelly et al. [2008] argue that, based on the exceptionally stable geochemistry of erupted lavas between 1972 and 2004, essentially constant magma temperatures must have existed without significant crystallization, magma intrusion and mixing. Indeed, they find that the mineral modes of recently erupted lava bombs are similar to those of the summit lavas erupted over the last ∼17 ka. Dunbar et al. [1994] emphasizes that crystallization has not been considered a significant process in the last 30 years. Based on different methods such as trace element diffusion modeling and ²³⁸U series geochronology, the most common age of the anorthoclase crystals is hundreds [e.g., Sumner, 2007] to thousands of years [Reagan et al., 1992]. Isotopic studies of Sr, Nd, Hf and Pb performed by Sims et al. [2008] show a relative compositional uniformity that also suggests that Erebus is in a steady state and that either its source is isotopically homogeneous, or the mixing of batches of parental basanite with different isotopic compositions is extremely efficient.

4.2. Idealized Magmatic System

[35] The aforementioned geological information guided us in framing an idealized magmatic system suitable for numerical simulations. Given the evidence that the Erebus shallow magma system has reached a steady state, we conducted simulations so as to achieve steady state. The geometry of this idealized system is based on a 2D cross-section referenced in Cartesian coordinates. This setup is a compromise between computationally very intensive 3D simulations and the forced symmetry of an axisymmetric domain. We fixed the lake diameter to 40 m and its depth to 20 m. We assume a constant heat flux between the magma and the vertical walls of the lake because snow and ice on the crater walls most likely perturb the thermal stability. This heat flux was calculated using a thermal conductivity estimated through a relationship from Zoth and Hänel [1988] (k_wall = 2.90 W/m °C) for a cold volcanic country rock of intermediate composition. In the absence of information on the density and heat capacity values of the matrix rock hosting the lava lake, we used common values for basalt, respectively ρ_wall = 2700 kg/m³ and C_pwall ∼ 840 J/kg °C, which yields κ_wall = k_wall/(ρ_wall C_pwall) ∼ 1.3 × 10⁻⁶ m²/s. The occasional Strombolian activity suggests that gas can be segregated at some level of the plumbing system. Although we do not simulate a gas phase, we use this information to derive additional knowledge about the conduit geometry. Jaupart and Vergniolle [1988] have argued that bubbles collect at the top of a chamber near the conduit and that collapse of the accumulated foam can generate gas slugs that rise through the conduit to burst at the surface [Chouet et al., 1997]. We thus simplified the plumbing system below the lake as a single vertical conduit linked to a magma reservoir. The location of seismic signals given by Aster et al. [2008] plus the depth of crystal formation [Dunbar et al., 1994], suggested fixing a conduit length of 400 m and a diameter from 4 to 10 m, within the range of widths reported in previous studies. Since the dimensions of the chamber are not constrained, we represented the magma reservoir as a 2-D domain in such a way that it is larger than the lava lake, 100 m wide and 50 m deep.

[36] The idealized 2D geometry (within the x-y plane) for Erebus is represented in Figure 4. The initial conditions over the whole computational domain are a uniform temperature of 1000°C and null velocity for all phases (melt, melt and crystals or mixture), and 30 vol.% of uniformly distributed anorthoclase crystals for bi-phase simulations. The pressure for the whole magmatic system is initially magma-static (mechanical pressure). The thermal boundary conditions (BC) were set for all phases so that the system is heated from below (BC1 = 1000°C) by a magma chamber, and cooled from above through the lava lake surface (BC8 = 0°C). All phases were kept at 1000°C on the reservoir sides (BC2 and BC3), and at 900°C along the conduit (BC4 and BC5). The vertical walls of the lava lake (BC6 and BC7) were set so that the heat flux through the walls would be of reduced intensity, C_m = −2°C/m. This gradient was obtained through equation (4) when the temperature gradient is established between a melt initially at 1000°C and the host rock at T_fwall = 0°C, which corresponds to cold country rocks in the far field (Figure 5). Even for these high initial temperature contrasts, the thermal evolution depends weakly on the initial contrast and is mostly governed by κ_wall and t. A series of simulations was performed in closed system, which means that no magma enters or exits the domain. In our closed system, we choose to fix the sidewall boundaries (BC2–7) as non-slip walls (NSW), while the top of the lava lake (BC8) and the bottom of the magma chamber (BC1) are free-slip walls (FSW). This last condition limits the effect that the arbitrary thickness has on the thermal evolution of the system. While this is a fair approximation for perfect mixtures of melt and crystals because of the absence of settling, the bottom boundary will cause crystal to accumulate for bi-phase simulations. The velocities for all phases approach zero on the sidewall boundaries, while the velocity gradients approach zero on the horizontal boundaries. One simulation was performed in open system, for which BCs are identical to those for the closed system, except that BC2 and BC8 were set as inlet and outlet boundaries, respectively (Figure 4). At the inlet, the bi-phase flow had a pressure of 11 MPa, a solid volumetric concentration of 30 vol.%, a phase-weighted horizontal speed of 2.3 × 10⁻⁵ m/s; at the outlet, the pressure was constant and equal to atmospheric (0.1 MPa).

[37] The size distribution and shape of crystals was also simplified. We calculated that a spherical crystal of 5 cm in diameter would have the same volume as a crystal of 6 × 3 × 3 cm. The spherical shape is a necessary assumption in order to compare the behavior of the bulk viscosity of the mixture with the bi-phase simulation in a consistent manner (see 2 and 2). We took anorthoclase properties such as density from the literature, and estimated other properties such as thermal conductivity and heat capacity based on the end-member composition of the anorthoclase crystals (albite). Table 3 summarizes all the physical properties for the crystals and the phonolite melt used in the simulations. The angle of internal friction between particles has, to our knowledge, never been determined in such viscous suspensions. After trying values of 0° and 24° without obtaining noticeable changes in average temperature and velocity patterns, we set this parameter to 0° so as to decrease calculation time. We chose an intermediate value for the particle restitution coefficient, e, of 0.8 (between cork, 0.6, and glass, 0.9). We neglected crystallization given the estimates of crystal ages in the literature and the limited evidence for its significance in the real system (see 11). Should it occur, the crystallization process would continuously grow existing crystals or nucleate new ones, which would eventually contribute to changes in the physical and kinetic properties of the magma (viscosity, density and drag). The validity of this assumption should be assessed in the future by comparing simulation results to the equilibrium conditions given by petrologic data for Erebus phonolite; such equilibrium conditions have not been published yet. We tested values for the minimum void fraction of a packed bed of crystals of 0.35 (as in Costa [2005]) and 0.37 [McGeary, 1961] without obtaining noticeable differences and we used a packed bed void fraction, ε_m*, of 0.35 as discussed later.

Table 3. Parameters Used in Physical Modeling of Erebus

Parameter	Unit	Crystals			Melt
		Symbol	Value	Reference	Symbol	Value	Reference
Density	kg/m3	ρs	2.59 × 103	Anthony et al. [1995]	ρo	2.455 × 103	Output from Conflowa
Thermal conductivity	W/m °C	ks	2.0	Clauser and Huenges [1995]	km	1.53	Clauser and Huenges [1995]
Heat Capacity	J/kg °C	Cps	2218.8	Navrosty [1995]	Cpm	1367	Output from Conflowa
Thermal diffusivity	m2/s	n/a	3.5 × 10−7b		κmb	4.56 × 10−7
Thermal expansion coefficient	1/°C	n/a	n/a		α	8.054 × 10−5	Output from Conflowa
Volume fraction	Dimensionless	εs	0.3	Kelly et al. [2008]	εm	0.7
Packed bed void fraction	dimensionless	n/a	n/a		εm*	0.35
Maximum packing fraction	dimensionless	n/a	0.65	Costa [2005]	n/a	n/a
Particle diameter	m	dp	5 × 10−2		n/a	n/a
Restitution coefficient between particles	dimensionless	e	0.8		n/a	n/a
Angle of internal friction	° (degree)	ϕ	0		n/a	n/a

• a Mastin and Ghiorso [2000].
• b Calculated value.

[38] We calculated the initial density ρ_o, the thermal expansion α and the heat capacity C_pm of the melt with Conflow [Mastin and Ghiorso, 2000], using the following parameters: a melt composition of anhydrous whole-rock of phonolite lavas given by Kelly et al. [2008], 0.26 wt.% of H₂O, a magma chamber pressure of 11 MPa, 30 vol.% of anorthoclase crystals, and a temperature of 1000°C. The parameters α, C_pm and k_m were kept constant (Table 3). The obtained value of ρ_o (2455 kg/m³) is close to the estimation of Eschenbacher [1998] for Erebus phonolite. The thermal conductivity of the phonolite, k_m, was estimated through a relationship of Zoth and Hänel [1988] using a volcanic rock of intermediate composition at T_o = 1000°C. Melt viscosity is of particular importance because it controls the rate of transport of matter and, thus, of energy. We used the Newtonian model given by Giordano et al. [2008], which predicts the non-Arrhenian temperature and compositional dependence of melt dynamic viscosity, μ_m, for naturally occurring silicate melts at atmospheric pressure. This viscosity appears in the deviatoric term of the stress tensor in equation (30). Based on the mean of the major element concentrations of matrix glass separated from lava bombs [Kelly et al., 2008] plus the water content of melt inclusions in anorthoclase (0.26 wt.%) measured by Seaman et al. [2006], we have calculated the coefficients B and C from the following relationship given by Giordano et al. [2008]:

Where A = −4.55 (coefficient independent of composition), B = 8642.3, C = 319.6 and T_m is the melt temperature in Kelvin. The equation (11) is valid over a range of viscosity values of 0.1 to 10¹⁴ Pa s and temperatures of 245–1580°C [Giordano et al., 2008]. The typical melt viscosity for hot and cold Erebus magma in our simulations is of the order of 3.3 × 10⁴ and 10⁸ Pa s, respectively.

[39] Our simulations include the lake surface; it is thus necessary to determine the temperature at which cooling quenches the magma into a crust. Physically, this would be the glass transition temperature, T_g, which is the temperature separating the liquid (relaxed) from the glassy (unrelaxed) state. We decided to clip the value of viscosity at the temperature at which T_g occurs, so that the melt shows a strong viscosity increase that approximates sub-solidus conditions while the absolute value of viscosity still ensures numerical stability. The clipping method underestimates the viscosity that a colder glass can reach. The validity of this assumption will be discussed in 25. Experimental and theoretical approaches suggest a range of values for T_g. Giordano et al. [2005] used differential scanning calorimetry and a fixed cooling/heating rate of 20°C/min to find an experimental peak-value for T_g of 670°C for a Teide phonolite. Russell and Giordano [2005] predicted a value of T_g for the end-member melt compositions diopside-anorthite-albite of 720°C, 857°C and 796°C, respectively. The thermodynamic model of MELTS [Ghiorso and Sack, 1995] yields an estimate of 600°C for the solidus temperature of our melt. The factors A, B and C from the viscosity model of equation (11) can also be used to estimate T_g [Dingwell et al., 1993; Giordano et al., 2008]:

This provides an estimate of 570°C, which corresponds to a viscosity of 10¹² Pas. This viscosity value corresponds to a relaxation timescale for macroscopic melt properties of about 15 min and to cooling rates of 10°C/min [Giordano et al., 2008]. These considerations restrict the probable range for T_g between 600 and 800°C. Temperatures around 600°C correspond to viscosities of ∼10¹¹ Pa s, which cause convergence failures due to a large increase of pressure and velocities for some cells in the top of the lake. We thus explored two clipping temperatures: 700°C, which is the lowest possible value that still ensures numerical stability, and 800°C, which approximates the value of 796°C given by Russell and Giordano [2005]. We first quantified the cooling rates of magma in the lake for simulations considering free-slip condition at the top of the lake (Figure 6a) and the respective viscosities as a function of temperature (Figure 6b). We found that cooling rates can be better compared to experimental values when the top of the lake is set to non slip (Figure 6c). We also obtained a horizontal average of melt velocity for all the cells below T_g (Figure 6d). When T_g = 700°C, melt velocity decreases from 3 × 10⁻⁵ m/s at a depth of 14 m below surface to 10⁻⁷ m/s just below the surface of the lava lake, because the velocities approach zero at the horizontal boundaries. In contrast, melt velocities remain almost constant above 10⁻⁵ m/s when T_g = 800°C. We note that changing T_g from 800°C to 700°C has a greater influence on the velocity profile than changing the slip conditions. When T_g = 800°C, a thick quenched central instability forms after 69.4 days, sinks toward the deeper part of the lake, preventing thermal exchange between lake and conduit to occur and driving the whole lake to drop below T_g (i.e., reaching a plateau average viscosity of 8 × 10⁶ Pa s; Figure 6b). When T_g =700°C, the whole lake remains above T_g (Figure 6b). Keeping in mind that T_g is a temperature interval related to the cooling/heating rate rather than a unique temperature value, we chose a clipping temperature of 700°C because it yields the lowest velocities at the surface of the lake (Figure 6d) and it is close to the value found by Giordano et al. [2005] at similar instantaneous cooling rate (Figure 6c). The corresponding viscosity, 10⁸ Pa s, is just below the experimental measures of 10⁹–10¹³ Pa s for the glass transition of silicate melts reported by Webb and Knoche [1996] at cooling rates of 5°C/min and T_g temperatures spanning the wide range 400–900°C.

4.4. Results Concerning How Melt and Crystals Are Treated in the Simulations

[44] In this section we present the results for the three selected simulations with 4–m-conduit diameter in closed system. Features of the convection in the whole system for these simulations are detailed in Figures 7a–7u and early stages of the convection in the lava lake are described in Figures 8a–8r. Details of specific numerical adjustments are detailed at the beginning of 14–16.

4.4.1. Pure Melt

[45] We tested numerical accuracy by keeping all the physical parameters of the pure melt configuration constant, as stated in Table 3, and varying the grid size between 50 × 50, 25 × 25, and 10 × 10 cm². Table 4 shows that the quality of the solution improves as the grid is refined, although the residuals of the norms of U_m and V_m do not evolve along the same trend. The two finer grids have the smoothest temporal evolution of both the time step history (Figures A3a–A3c) and the temperature decay (Figures A3d–A3f). The solutions with the coarsest grid show only a few instabilities surrounding a large central instability dropping from the surface of the lava lake (Figure A3a). In the smallest grid, on the other hand, there are 17 thermal plumes originating from the lake surface (Figure A3c) and 50 plumes from the top of the magma chamber (Figure A3f). The run with 10 × 10 cm² cells thus provides a better resolution of the RT instabilities and was consequently selected. It corresponds to 1000 × 4700 computational cells. Numerical accuracy was checked by varying the tolerance from 1 to 0.5 and 0.1. The solutions yield very similar evolutions of the temperature profiles, which indicate that the largest value of tolerance is sufficient.

[46] In 6, we showed that the number of instabilities varies inversely with viscosity for a simulation with given BCs and initial conditions. We used Odé [1966] and Lister and Kerr [1989] to predict λ_c (equation (9)) and consequently the number of instabilities for those simulations. These theories predict that 17 and 52 instabilities should form at thermal boundary layers with thicknesses h_c of 0.8 and 0.65 m for the lake and the chamber, respectively. These values are comparable to the results of our simulations, which lead to 17 and 50 instabilities for lake and chamber, respectively. The presence of the conduit favors the merging of the instabilities at the junction between conduit and chamber. It also opens a window in the ceiling of the chamber where instabilities would otherwise have formed. This layout explains the 4% discrepancy between the theoretical values and those from our simulations.

[47] We turn next to the downward motion of the instabilities. For a melt with ΔT = 900°C and Ra = 9.5 × 10⁹ (lake) or ΔT = 100°C and Ra =1.7 × 10¹⁰ (chamber), temperature-driven convection typically yields Rayleigh-Taylor instabilities [Marsh, 1988], observed as a progressive front when the downgoing cold plumes entrain more and more hot fluid with time (Figures 7a–7f and 8a–8e). Less than 15 days after the beginning of the simulation, we observe the formation of a crust atop of the lake, which is marked by the clipped viscosity that represents the glassy state. There are no rising plumes, but, with time, the cold and hot fluids inside the lake and chamber eventually mix. The RT instabilities take 4.8 years to reach the bottom of the lake, while the plumes forming at the top of the chamber do not reach its base within the 30 years of the simulation (Figures 7c and 7f). After the first year, the flow pattern quickly reaches a steady average temperature in the magma chamber. In the lava lake, this process slows continuously to reach a homogeneous temperature of ∼950°C around year 10 (Figure 7g). The flow pattern in the conduit always shows a very gradual downward motion (Figures 7d–7f) because 30 years is not long enough to induce a significant return flow that can offset the cooling effect of the 900°C conduit walls.

4.4.2. Crystals as a Part of the Melt

[48] Based on the vertical profile of the average temperature, the simulations with different grid sizes show only slight differences in the lava lake before year 10 of the simulation; the simulation obtained with the finest grid being warmer by <10%. After 10 years, the three grid sizes converge on the same temperature range (around 500°C), which means that the lava lake completely freezes. No temperature differences are found in the magma chamber and they amount to 2% in the conduit. We selected the simulation with a grid size of 25 × 25 cm² and changed the tolerance from 1.0 to 0.5. Residuals are of the same order (Table 4) and vertical temperature profiles are similar, suggesting that the largest value of tolerance is sufficient. This simulation shows a cold crust of 1.75 m (Figure 8f); using this value in equation (9) yields an estimate λ_c and indicates the formation of 8 instabilities, consistent with the number of instabilities we observe (Figure 8f).

[49] The cold crust at the top of the lake develops a central instability that sinks (Figures 8f–8h) and reaches the shallower part of the conduit after 3.5 months (Figure 8i). Until then, the convection pattern is characterized by an ascending core of fluid through the central part of the conduit and descending fluid at the sides. The extension of this flow pattern to a 3D geometry would be a core-annular flow, a terminology we use hereafter for convenience. This type of convection ceases when the central instability starts folding on itself like a periodic disturbance to create a “plug” that seals the lake at month 6. This plug (Figures 7h and 8j) leads to a continuous cooling of the lake (Figures 7i and 7j). All the heat in the system is then lost by constant cooling of the conduit, which features continuous downward motion and two independent convecting systems develop in the lake and the chamber (Figures 7k–7m). The chamber has the highest convective rates with a velocity oscillating around 7 × 10⁻⁵ m/s. The temperature of the magma chamber and the conduit remain constant around 990°C and 900°C, respectively (Figure 7n). Most of the lake reaches glass transition temperature after the first year of the simulation and after year 5 the lake reaches a temperature of about 500°C.

4.4.3. Crystals as a Separate Phase

[50] The only difference between grid sizes of 50 × 50 and 25 × 25 cm² is the definition of the RT instabilities at the onset of convection. With the 25 × 25 cm² grid, we counted at least five instabilities (including the central one) that are only poorly resolved by the coarser grid. Two lateral instabilities progress along the vertical walls of the lake to reach its bottom with a small delay compared to the central instability. These lateral instabilities are progressively mixed within the lake during the first year of simulation, which is characterized by highly unsteady convection at all grid size. These instabilities do appear in the grid size (50 × 50 cm²) that we selected to conduct the tolerance study (Table 4). This choice was motivated by the fact that these simulations are time-consuming; one simulation with the coarse grid running on 144 processors takes 15 days to complete 30 years of simulated time. We found that the runs are almost insensitive to the tolerance values because the corresponding average temperatures differ by only ∼2%. The time step history is the smoothest when tolerance equals 1.0, so we selected that run even if the residuals given by the other runs were slightly smaller (Table 4). This run shows a crust of 2.5 m (Figure 8k), from which equation (9) allows us to predict 5 instabilities. This is the number of instabilities we observe (Figure 8k).

[51] As instabilities sink deeper in the lake, they merge with each other and smear out (Figures 8k and 8l), except the central instability, which sinks toward the shallower part of the conduit (Figures 8l–8n) that it reaches after 2.5 months to form a periodic folding disturbance, the natural frequency of which can be compared with analogue models (Figure A4 and details of this in 25). The periodic folding instability behaves differently to that observed in the simulation with crystals as part of the melt. The folding partly plugs the upward flow until year 1 before being displaced laterally by the continuous flow coming from the conduit (Figure 8o). After year 1, the fluid in the lake keeps convecting and partly chokes the conduit (Figures 7o–7q), and the convective rate decreases. The flow pattern in the conduit evolves from a core annular flow (Figure 7r) into vertically stratified flow (Figure 7s) and later into upwards-only motion (Figure 7t). The continuous cooling of the system isolates the lake and part of the conduit from the rest of the system, which then features two sustained core annular flows from chamber to the lower part of the conduit and from lake to the shallower part of the conduit, respectively (Figures 7s and 7t). As in the simulation with crystals as part of the melt, most of the lake reaches glass transition temperature after the first year of the simulation (Figure 7o) and attains a temperature of about 500°C after year 5 (Figure 7u).

4.4.4. Comparison of the Three Scenarios for Representing Phases

[52] A striking difference between the three configurations is the convection inside the conduit and the flow velocities in the lake. The bi-phase configuration shows a core-annular and a vertically stratified magma flow inside the conduit that facilitates thermal exchange between the lake and the magma chamber. When crystals are considered as a part of the melt, a similar series of patterns occurs, with an earlier shift between a single and multiple core-annular flows. In contrast, the pure melt simulation shows melt being pushed down the conduit by the front of RT instabilities descending within the lake, without generating a bi-directional flow. We surmise that a return flow would be established at much longer times. Within the time span of our simulations, however, the presence of core-annular flow along the conduit is enhanced by the presence of crystals.

[53] The velocity of the fluid is higher in the lower section of the lake during the first month (Figure 8p), and the fast-moving front progressively migrates upwards (Figures 8q and 8r) under the influence of both thermal convection and crystal settling. The presence of crystals influences the convective rate, despite the fact that particle movement is well coupled to the liquid (Figures 8p–8r). When crystals are present as a separate phase, the initial convective rate over the first 2 months is relatively high (the highest of the three) (compare Figures 8p–8r). The central instability temporally plugs the conduit at 2.5 months (Figure A4), which temporally isolates the lava lake from the conduit and decreases the convective rate (Figure 8r). In a mixture with crystals as part of the melt, it takes slightly longer (3.5 months) for the instability to reach and plug the top of the conduit (compare Figure 8i with Figure 8m). The pure melt behaves differently; after 1 month only a weak change in velocity occurs at the base of the lake (Figure 8p). The instabilities descend progressively from the top of the lake to reach only a few meters by month 3 (Figures 8a–8c and 8p–8r). The presence of crystals entrained by the central instability enhances the velocity at which the bi-phase flow sinks at the vertical of the conduit (compare Figure 8m with Figure 8h). The presence of crystals and the way they are associated with the melt have a lasting influence on the flow velocity at the surface of the lake, as discussed further in 23.

4.5. The Influence of Conduit Dimensions

[54] We examine the influence of the conduit diameter while keeping the other parameters of the three different configurations constant. Figures 9a–9u show temperature and vertical velocities for the “pure melt,” “crystals as a part of the melt,” and “crystals as a separate phase” configurations with a 10-m diameter conduit. The comparison between Figures 7a–7g and 9a–9g suggests that increasing the conduit diameter from 4 to 10 m for a pure melt does not affect temperature nor flow velocities in the lake and chamber. The temperature remains higher in a 10 m diameter conduit, where a core annular flow develops unlike in a 4 m conduit, where we observe a downward motion. For a 4-m-diameter conduit, crystals as a separate phase sustains the heat transfer from chamber to the lava lake for the longest period (compare Figures 7k and 7l with Figures 7r and 7s). Increasing the conduit diameter increases this period even further. The overall temperature of the system increases and remains above glass transition (compare Figure 7u with Figure 9u). The core-annular flow is maintained for all the configurations and all simulated times (compare Figures 7a–7u with Figures 9a–9u, respectively). This is mostly due to the inability of the central instability to plug such a wide conduit. A 10-m diameter conduit thus allows the full system to approach steady state convection.

4.6. Sensitivity to Boundary Conditions

[55] The numerical scheme of MFIX is based on finite volumes. When the melt cools and contracts, the imposed closed system causes mass loss. This loss is small over the 30 years of simulations (∼1%). The algorithm of MFIX, however, always enforces equation (14) at the end of each iteration, which causes the small mass loss to be translated into a crystal volume fraction loss. In our bi-phase simulations, the densities of crystals and melt are close to each other. As a result, the crystal volume fraction loss becomes non negligible. We estimate that it is equal to 5% in the chamber and conduit. Up to 8 vol.% crystals are thus lost in the lake because the temperature change is higher there than in the rest of the system. Simulations in open system are not subject to these mass and crystal load losses. Whether or not these losses affect convection dynamics can be assessed by comparing the bi-phase simulation with 10-m conduit diameter and in closed system with a similar simulation but in open system. We first provide a mechanical description of the closed system simulation with crystals as a separate phase and a 10-m conduit diameter before highlighting the differences introduced by the open system. Figures 10a–10c, 11a–11t, and 12a–12g highlight the effect of changing the boundaries of the system from closed to open. Figure 10 shows the evolution of temperature and bulk densities in the three regions of our closed system; the main stages of the convection and crystal settling within the lake for a 10-m-diameter conduit system are shown in Figures 11a–11e for a closed system and in Figures 11f–11j for an open system. Similarly, we compare the stages of convection for the chamber in a closed system in Figures 11k–11o to those for an open system in Figures 11p–11t. Details of the evolution of temperature variations in the three regions and crystal content inside the conduit for an open system are shown in Figure 12.

4.6.1. Crystal Settling in a Closed System (No Feeding)

[56] Over the 30 year period, the magma cools down by 266, 89 and 60°C in the lake, conduit, and chamber, respectively (Figure 10a), and one can expect its bulk density to increase. Our simulation indicates, however, a decrease of bulk density of the magma in the lake, conduit and the chamber (Figure 10b). This is a likely result of crystal settling in the system. The steady state values of the bulk density, from ∼4 to 30 years, correspond to the cooling of the three domains. These results allow us to describe mechanistically the main stages of the magmatic system evolution.

[57] At the beginning of the simulation, the temperature of the vicinity of the surface of the lake decreases (Figure 10a), resulting in the formation of a thin crust of ∼1.5 m and in the growth of the first instability at the center of lake after ∼15 days (tag 1 in Figure 10c). Instabilities carry a few crystals enveloped by the quenched liquid (∼20 vol.%) into the shallower part of the conduit (tag 2 in Figure 10c). Both crystals and the melt reach a maximum speed of 5 × 10⁻⁵ m/s. This strong cooling causes the bulk density of the lake to rise rapidly over the first 3 months (Figure 10b). In the conduit and the chamber, the bulk density does not increase as much as in the lake because these two domains are insulated. Unlike the 4-m conduit simulation, the plugging of the conduit-lake junction does not occur (tag 3 in Figure 10c). At month 3, the central plume shifts laterally and merges at month 6 with one of the lateral instabilities that has been simultaneously developing. This partially isolates the lake (tag 4 in Figure 10c), causing a convective pattern that keeps the volume of crystals approximately constant (∼17 vol.%) inside the lake (Figures 11b and 11c). A hotter plume rises from the conduit at ∼1.5 yr and progressively destroys the partial isolation of the lake by generating stronger convection in the lake (tag 5 in Figures 10a and 10c). This prevents the system from cooling further but sedimentation continues until year ∼10 to stabilize at 6 vol.% of crystals in suspension (Figures 11d, 11e, 11n, and 11o). This temperature increase marks the progressive transition toward steady state at year ∼4 (Figure 10a), at which point the center of the lake remains at a temperature above glass transition (734°C).

[58] The crystal content in the magma chamber never becomes homogenous during the first unsteady 4 years because cool and crystal-poor liquid keeps pouring from the conduit into the chamber. A crystal-rich (60%) layer appears at the chamber bottom when steady state is reached, but it cannot be described as “stagnant” because this layer is eventually disturbed by the convection generated by the arrival of cooler liquid from the conduit, containing less crystal. Remains of this layer can be seen in Figures 11n and 11o. The bottom layer remains enriched in crystals while the top of the magma chamber becomes almost crystal-free (Figure 11o). Some cells briefly reach ε_m* or 1% above because of plasticity, regardless of the value of ε_m* (35 vol.% or 37 vol.%).

[59] Taking the system as a whole, the convective pattern in steady state remains stable until the end of our 30 year simulation. Relatively crystal-rich, hot fluid reaches the lake from the conduit, sustaining the convection pattern and 6 vol.% of crystals remain suspended within the convecting magma. Should all the crystals that correspond to the 24 vol.% that are settling from the three regions accumulate at the bottom of the chamber and remain there at maximum packing fraction (1 − ε_m* = 65 vol.%), the crystal-rich layer would reach 11 m in thickness. Two large and steady convective cells are observed in the lake (Figures 11d and 11e) while chamber convection does not show any characteristic pattern over time except at the end of the simulation, when two regular and symmetrical cells are established, the one on the right side circulating in an anticlockwise direction the other circulating clockwise (Figure 11o).

4.6.2. Crystal Settling in an Open System (Permanent Feeding)

[60] Based on changes in temperature profiles we differentiate two main states: unsteady behavior in which a maximum temperature drop reaches a plateau before rising again to reach steady state by year ∼4. In steady state, the open-system temperatures drops by 79, 82 and 46°C in the lake, conduit, and chamber, respectively (Figure 12a). This is warmer (compare Figure 12a with Figure 10a) and more crystal-rich than the closed system (compare Figures 11f–11j with Figures 11a–11e for the lake, and Figures 11p–11t with Figures 11k–11o, for the chamber). Temperatures in both the shallow and the deep parts of the conduit decrease progressively during early stages before increasing at later times (Figures 12a–12d). This translates into a small variation of the shallow temperature contrast, from 17 to 8 to 14°C, while the contrast of deeper-seated temperatures has larger variations, from 56 to 35 to 45°C. Details of the evolution of the crystal content inside the conduit are shown in Figures 12e–12g. Ultimately, shallow and deep temperatures in the conduit increase as a consequence of a fully developed convection between the deep and hot region and the shallow and cooler region.

[61] In the chamber, several regular patterns of rolling convective cells appear at an early stage of the open system simulation (Figure 11p), before resolving into fewer cells as the mixture becomes more homogenous (Figures 11q and 11r) and as a crystal-rich layer appears at the chamber base (Figures 11s and 11t). Crystal-melt segregation at the base of the chamber is enhanced as crystals deposit, and manifests itself in the form of tiny plumes rising from a compacted region where hindered settling occurs (Figure 11p). Fluid motion generated by the convection in the chamber is practically horizontal right above the compacted region at later stages of the simulation (Figure 11t). This is in contrast to the closed system, where the much slower convection barely shows any characteristic pattern over time except at the very end of the simulation (compare Figures 11k–11n with Figures 11p–11s, respectively). Another difference with the closed system is the 10-m-thick, crystal-rich (∼33 vol.%) layer that remains at the chamber bottom throughout the steady state (compare Figures 11n and 11o with Figures 11s and 11t, respectively). The cold flow pouring from the conduit sustains convection in the upper part of the chamber without disturbing this compact layer, the temperature of which remains slightly higher than the rest of the domain (Figures 12c and 12d) and causes enhanced melt extraction (Figure 11s). The thickness of accumulated crystals remains around 10 m whereas the respective fractions of melt and crystals vary due to re-entrainment of crystals in the chamber (Figure 11t). The arbitrary thickness of the chamber conditions the thickness of the settling horizon and the convective patterns of the chamber. Both closed and open systems display a crystal-poor layer at the top of the magma chamber on each side of the entrance to the conduit (compare Figures 11n and 11o with Figures 11s and 11t). This may suggest that the effects of the thickness of the chamber are probably limited to this level of the reservoir.

[62] In open system, a crystal-rich layer forms at the bottom of the lake on each side of the conduit (Figures 11i and 11j), mirroring that of the chamber base. Both closed and open systems result, however, in similar convection pattern in the lake. Motion is clockwise at the left side of the cell and counterclockwise on the right side of the cell (Figures 11b, 11c, and 11h), a pattern that will be reversed over time once steady state has been reached (e.g., Figures 11d, 11e, 11i, and 11j).

[63] Figure 13 describes characteristic flow patterns occurring in a 10-m conduit inside a 10-m-diameter conduit, which we call respectively sinusoidal (Figure 13a), core annular (Figure 13b) and vertical stratified flows (Figure 13c). The convective flow in a closed system presents a sinusoidal profile combined with a vertically stratified flow throughout the entire period of simulation (Figures 9r–9t and 13c). Open system causes changes in flow patterns; after several months, the conduit develops a core annular flow (Figure 11h) that evolves into a sinusoidal pattern at ∼year 1 (Figures 13a). Later in the simulation (e.g., year 30), the open-system conduit features a well-established core-annular flow (Figures 11i and 11j), replacing the sinusoidal flow. We observe at year 2.6 that the sinusoidal flow is caused by the interaction of cold and hot fluids. The cold fluid (903°C) transports a magma with relatively higher crystal content (23 vol.%) and moves downward along the walls of the conduit. The hot fluid (911°C) transports a lower crystal content (19 vol.%) and ascends in the core of the conduit, from the chamber to the surface of the lake (Figures 12c and 12f). The sinusoidal regime displays a longer and longer wavelength as this instability approaches the top of the chamber (Figure 12c). It can be observed in the deeper and shallower part of the conduit that the contrast between the temperatures of ascending and descending fluids decreases to reach its lowest value at year 2.6 (Figure 12c). Then the temperature of the ascending fluid rises again (Figure 12d). The difference of crystal contents between ascending and descending fluids remains constant within 1 vol.% in the lower part of the conduit while a maximum difference of up to 4 vol.% is observed in the upper part of the conduit (Figures 12e–12g); the absolute value of crystal content increases slightly (2 vol.%) between year 2.6 and 26 in the lower part of the conduit. The maximum observed contrast in temperature is 56°C (Figure 12b).

4.7. Analysis of Salient Model Outputs: Behavior of the Lake Surface

[64] We compared average surface velocities over a depth of 1 m for all seven simulations (Figure 14a). We also used the method of Welch [1967] to compute the Power Spectral Density (PSD) of the velocity average time series to compare the main frequency peaks. This estimate was achieved by removing the linear trend and estimating the PSD with 12 overlapping windows (50% overlap between segments) to improve the signal-to-noise ratio of the spectrum. The spectral amplitude was normalized and the main frequencies shown in Figure 14b are below Nyquist frequency, which is equal to 4 × 10⁻⁷ Hz for our sampling rate of 15 days.

[65] Overall, the trend of the average velocity is identical for both conduit sizes and the three different configurations tested to study the Erebus (Figure 14). For a closed system, the magnitude of the velocity depends more on whether the crystals are considered as part of the melt or as separate phase than on conduit diameter. The velocities of the crystals as part of the melt fall between the velocities of the bi-phase and pure melt simulations. Among all simulations, pure melt magma always shows the lower surface velocities, with an a average value of 8.6 × 10⁻⁸ m/s early on in the simulations (over the first 3 months) to reach a constant value of a 6.6 × 10⁻⁸ m/s over the rest of the period. This behavior is similar for both conduit diameters. The velocity of magma with crystals as part of the melt, in a closed system, with a 10-m conduit increases from 3.2 × 10⁻⁷ (day 15) to 3.3 × 10⁻⁶ m/s during the first half-year before dropping continuously until about year 11, when it starts to oscillate around 9.8 × 10⁻⁷ m/s. When the conduit diameter is only 4 m, this behavior is similar until month 3.5 when the velocity reaches a peak of 1.8 × 10⁻⁶ m/s before decreasing due to plugging of the conduit. It then increases again slightly (up to 1.5 yr) to reach a constant value similar to that in a 10-m conduit but without presenting any oscillating pattern. Liquid carrying crystals as a separate phase, in a closed system, with a 10-m conduit shows an increasing trend of velocity during ∼1.4 yr, from 2 × 10⁻⁷ (day 15) to 4.9 × 10⁻⁶ m/s. The surface velocity reaches later an average constant value of 3.2 × 10⁻⁶ m/s. When the conduit diameter is only 4 m, the velocity behavior is similar over the first year, when it decreases from 6 × 10⁻⁶ to 2.7 × 10⁻⁶ m/s over the next 17 years. It then increases sharply to 3.7 × 10⁻⁶ m/s and remains at that value for the rest of the simulated period. In an open system with crystals as a separate phase and a 10-m conduit, the surface velocity is higher than in closed system. It sharply decreases from 1.7 × 10⁻⁴ (day 15) to 2.2 × 10⁻⁵ m/s at year 4, a value that remains approximately constant until year 21 when it increases smoothly to an average value of 6.9 × 10⁻⁵ m/s (from year 26 to 30).

[66] Whereas a pure melt features a single large initial amplitude of 10⁻⁷ m/s (0–3 months), the other configurations display velocity oscillations. The velocity of the mixture oscillates around ∼10⁻⁶ m/s with a period of 0.4 yr (over both the entire time of the simulation and from 11 to 30 yr) and that of the bi-phase in a closed system around 3.2 × 10⁻⁶ m/s with a period of ∼0.3 yr (over both the entire time of the simulation and from 1 to 30 yr). For the bi-phase simulation in an open system, the surface velocity oscillates around ∼10⁻⁵ m/s with a period of 0.1 yr from 4 to 21 yr and over the last 4 years. Therefore, the velocity of the lake fed by 10-m-diameter conduit increases and decreases every half a period, or every 3, 2, and 1 months for a mixture in closed system, a bi-phase magma in closed system, and a bi-phase magma in open system, respectively (Figure 14b).

5. Discussion

[67] We carried out seven numerical simulations of an idealized magmatic system broadly aimed at representing the Erebus lava lake. Four simulations were single-phase and three were bi-phase, which allowed us to test the effects of crystals on convection while changing conduit size from 4 to 10 m and the system boundaries from closed to open. Our bi-phase simulations with a sufficiently large conduit (10 m) show that convective motions under steady state conditions are sufficiently vigorous to keep a proportion of the crystals in suspension and well-mixed. In particular, when the closed system reaches steady state around year 4, 65% of the initial 30 vol.% crystals have settled. The dense layer formed at the bottom of the chamber is then disturbed by convection. Complete settling is thus never achieved in any region of the system, leaving 6 vol.% crystals in suspension from year 15 until 30 and presumably for much longer. In the open system, only 34% of the initial crystal content settles, leaving about 20 vol.% in suspension in late steady state. Simulating magma by a mixture with crystals as part of the melt, on the other hand, has the major issue that crystal settling is not taken in account. Crystals, if considered as a separate phase, play a fundamental role in increasing apparent viscosity locally and creating complex circulation pattern, and in enhancing the vigor of fluid dynamics. However, the mixture theory can provide some features of the system, such as (i) the formation of a central instability, (ii) the average temperature evolution, and (iii) the average velocity range of the flow motion at the surface of the lake.

[68] Crystal size and the contrast between crystal and melt densities (Δρ = ρ_s − ρ_m) are the two main controls on the amount of crystals remaining in suspension in the convecting magma. Using Stokes' velocity, v_s, we calculated settling velocities of 10⁻⁶ m/s for the closed system, a value of the same order of magnitude as the flow velocities given by our bi-phase simulations. This balance between settling and convection is in qualitative agreement with the theoretical findings of Huppert and Sparks [1988] and Burgisser et al. [2005]. Our simulations also show that some crystals come out of suspension and accumulate on the lake and chamber floors, as reported in other studies of crystal settling [Martin and Nokes, 1988, 1989; Martin, 1990; Worster et al., 1990].

[69] Lavorel and Le Bars [2009] studied the sedimentation of spheres in a convective medium following the work of Martin and Nokes [1988, 1989]. They observed that the solid fraction remaining in suspension is a function of Ra and of the density ratio between particles and the fluid (Δρ/ρ_m). Consequently, the number of particles N decreases as a function of time t as follows:

Where is the initial number of particles per cell (=38 particles/cell), t_conv is the time necessary to establish the temperature contrast inside the boundary layer, hence t_conv = (H²/π κ_m)(2Ra/Ra_c)^2/3 ≈ 10⁴s. N_eq is a constant equilibrium value related to the number of particles when the statistical steady state is reached (dN/dt = 0), v_s is the Stokes velocity and H is the height of the system (470 m). The scaling law characterizing that equilibrium constant is written as N_eq = σ(18π²)(κ_m²ν_aν/H⁴g²d_p²)(Ra/Ra_c)^4/3(Δρ/ρ_m)⁻²(S_cell/S_crystal), where σ is a constant or “efficiency factor” of order 1, κ_m is the thermal diffusivity of the magma, ν_a is the apparent viscosity which depends on the kinematic viscosity of the mixture ν as ν_a ∼ (Ra/Ra_c)^0.06ν, g is the gravitational constant, d_p is the crystal diameter (5 × 10⁻² m), S_cell (0.25 m²) and S_crystal (2 × 10⁻³ m²) are the surfaces of the cell and crystal, respectively. Thus, for Ra/Ra_c = 6 × 10⁹ and ν = 35.81 m²/s, ν_a = 137.8 m²/s. We applied equation (13) for a constant temperature (therefore at constant density and viscosity) and compared the results with those for our bi-phase simulation with a 10-m conduit (Figure 15a). Our data in closed system show that, after a sharp early decrease, the number of particles in suspension stabilizes around year 15 to reach a constant value of 8 particles per cell (∼6 vol.%). This profile matches the theoretical prediction over the first year, after which the number of particles in suspension keeps decreasing slightly slower than for the theory where the plateau corresponds to 7 crystals per cell. This difference is partly explained by the fact that the experiments of Lavorel and Le Bars [2009] were conducted in a rectangular tank with set non-slip wall conditions, whereas we have set different boundary conditions. Another factor is that our geometry is more complex than a single tank, which causes settling in the chamber and the lake to be affected by conduit circulation.

[70] Unlike the experiments of Lavorel and Le Bars [2009], the open system simulation has a permanent inlet/outlet of mass that limits comparison between our results and their theory. Our open-system simulation shows that most particles settle during the first 4 years to reach a plateau of 18 particles per cell from year 4 to 16 (Figure 15a). At the steady state (Figure 12a), a stepwise increase of suspended particles is noticeable after year 15 (Figure 15a). This is because a surge of temperature in the chamber around year 15 results in a delayed increase of the temperature in both the conduit and the lake until year 26, from which time it remains stable until year 30 (Figure 12a). This change in temperature is associated with a re-entrainment of suspended particles, the concentration of which increases until it reaches a new plateau at year 26 (Figure 15a). These variations, in turn, cause a velocity increase at the surface of the lake (Figure 14a).

[71] There is another way to understand crystal settling, which is mostly apparent in the open system simulation. Following the numerical study of crystals settling in convecting magma chambers by Verhoeven and Schmalzl [2009], and focusing on changes in temperature profiles and in particle concentrations in the chamber (Figures 15b and 15c), we differentiate two regimes of convection that affect the whole system (Figure 15a). These periods correspond to the unsteady and steady states that we determined for this simulation (Figure 12a). The thermal-dominated regime (T-regime) prevails up to year 4, a period over which crystals are mostly advected by the melt because the thermal driving force is strong enough to compete with gravitational settling forces. In this regime, the temperature and crystal content of the system decrease (Figures 12a and 15a) until particles have settled to form a ∼10 m crystal-rich layer at the base of the chamber (Figures 15b and 15c). This regime then evolves into a settling-dominated regime (C-regime, from year 4 to 30), which is characterized by particle-driven convection that keeps the chamber segregated into two layers. The bottom layer is warm and mainly consists of crystals, whereas the top layer is colder and has fewer particles in suspension (Figures 15b and 15c). In the C-regime, most crystals are permanently held in suspension in the whole domain (Figure 15a). The thickness of the crystal-rich layer remains constant while the temperature sometimes displays a negative gradient across this layer (Figures 15b and 15c). This reflects the association of an increase of temperature in the crystal-rich region with an increase in crystal concentration. Such changes are similar to those described by Verhoeven and Schmalzl [2009].

[72] The bi-phase simulations show intricate convective patterns. Patterns change as short-lived, multicellular convection in the whole system during the initial unsteady phase that later merge progressively into two large cells, first in the lake after steady state has been reached (closed and open systems, Figures 11d, 11e, 11i, and 11j), and then in the chamber at a later stage (closed system, Figure 11o). A similar phenomenon is seen along the conduit where convective cells increase in size while decreasing in number. Crystals tend to remain with a given cell but their thermal history is likely to reflect rather complex trajectories when the number of convective cells changes. The richness of convective patterns produced by our simulations has counterparts in the literature. Our bi-phase simulations in either open or closed system reproduce the particle movement photographed by Weinstein et al. [1988] when steady state has been reached in the lava lake. In both studies, particle paths (streamlines) are nearly parallel at the center of the lake, traveling fastest when the flow is mostly vertical. Weinstein et al. [1988] have shown experimentally that, even after most particles have settled, some of them could be indefinitely retained in the laminar convecting viscous flow. Likewise, our simulations reach a steady crystal load once a balance between the carrying capacity of the convection currents and settling has been established. Dufek and Bachmann [2010] observed in numerical simulations that crystals accumulating on the floor of a crystallizing magmatic reservoir form a layer where melt extraction is most efficient. Unlike Dufek and Bachmann [2010], our open-system simulation does not take crystallization into account but shows a similar behavior, whereby small fluid-rich plumes rise through a crystal-rich layer at early times (Figure 11p) while thermal convection re-entrains crystals having reached the top of this layer (Figure 11t). In the conduit, our simulations generally show that convection is characterized by an annular or vertically stratified flow, as experimentally reproduced by Beckett et al. [2011]. In the early stages of the convection in open system and during most of the closed system simulation, our results show a sinusoidal pattern due to the contact between upwelling hot and less crystal-rich magma confined to the central axis of conduit and annular descending cold, denser and crystal-rich magma that partially adheres to the sidewall of the conduit. This pattern has been experimentally shown by Huppert and Hallworth [2007]. The convective interaction of the annular current with the magma chamber, leaking down from the lake through the conduit and returning from the chamber as a central ascending flow, was also observed in experiments with a similar geometry (two reservoirs connected by a duct) [Bazarov et al., 2007]. This similarity explains why the structure of our convection currents is in good agreement with those described by Bazarov et al. [2007] for the overall magmatic system, except for the open-system simulation, which has a more complicated convective pattern.

[73] For a 4-m-diameter conduit and crystals as a separate phase, we observed that the tip of the central plume folds as it reaches the entrance of the conduit. Considering that the central plume has temperatures below that of glass transition and could therefore behave like a brittle solid, is this behavior realistic? Our assumption is that melt below 700°C behaves with a high (10⁸ Pa s) but finite viscosity. Brittle behavior appears in a melt with such viscosity if it is strained faster than ∼2 s⁻¹ (= ϑG/μ_m, with ϑ = 10⁻² as experimental constant and shear modulus G = 25 GPa [Webb and Dingwell, 1990]). The glassy center of the instability is sheared at typical rates ≤10⁻⁶ s⁻¹ (Figure 6d), which means that instability viscous folding is consistent with the capping of viscosity values. Glass transition occurs when strain rates exceed ∼10⁻⁶ s⁻¹ for melt viscosities of 10¹⁴ Pa s, which correspond to 245°C [Giordano et al., 2008]. Since such low temperatures are unlikely to be preserved very deep within a convecting lake, these considerations suggest that, to first order, our treatment of cold instabilities within the lake is realistic. Loubet et al. [2009] experimentally studied how a thin sheet of highly viscous corn syrup deforms when it falls vertically from a distance H through less viscous syrup onto an impermeable surface. They observed that, for a sheet/ambient viscosity ratio larger than 10 and ΔT > 300°, the sheet undergoes a folding instability. For a calculated ΔT of 475° and a viscosity ratio around 10⁴, MFIX reproduces this folding instability remarkably well. The corresponding scaling law for the folded sheet amplitude is δ = 0.41(‖v‖_m/f) where f is the frequency at which the instability oscillates (f = ‖v‖_m/H), and ‖v‖_m is the velocity of the flow in the bent zone. Our simulation yields ‖v‖_m = 2 × 10⁻⁵ m/s for a lake depth of 20 m. Therefore f = 10⁻⁶ cycles/s, which yields an instability amplitude δ of 8 m. This is the value we measured in the simulation (Figure A4).

[74] As expected, the result of some of our most restrictive assumptions, such as neglecting gas bubbles or crystallization, is to limit the range of agreement between our simulations and the natural system of Erebus. One of our motivations was to assess if thermal convection alone (without the added buoyancy of a gas phase) could prevent the lava lake from freezing. We found that a pure melt system would not freeze over 30 yr, regardless of conduit diameter. A lake containing crystals - be they treated as a separate phase or not - would see its average temperature reach the glass transition after 1 year if it is connected to the chamber by a narrow, 4-m-diameter conduit. A 10-m-diameter conduit, however, ensures that the average lake temperature remains above glass transition. Our results can be extrapolated to the natural systems qualitatively; they suggest that, given a wide enough conduit, the lake does not freeze and maintains live connection in presence of crystals. Another relevant result is that our 10-m conduit simulation in open system yields rather homogenous temperatures within the system. Such homogeneity is also observed at Erebus because the matrix glass composition, magmatic temperature (∼1000°C) and oxygen fugacity are time-invariant in the lake based on geochemical analyses of matrix glass [Kelly et al., 2008]. In the conduit of this simulation, the maximum temperature difference (56°C) between upwelling and downwelling magmas during the unsteady state, when there is still 30 vol.% crystals inside the conduit, is close to the value (65°C) reported by Sweeney et al. [2008] for the same crystal content. Our range of bulk viscosities (10⁴ to 10⁸ Pa s) for the shallower part of the magma system is consistent with the values reported by Sweeney et al. [2008] and the thickness (1.5 m) of the crust at the top of the lake given by our closed system is consistent with visual observations of a thin crust [Oppenheimer et al., 2009].

[75] The velocities of the magma at the surface of the lava lake calculated by our bi-phase simulations with a 10-m-diameter conduit are on the order of 10⁻⁶ m/s in closed system and 10⁻⁵ m/s in open system, which is 10⁵ and 10⁴ times smaller, respectively, than surface speeds of the Erebus lava lake estimated by Oppenheimer et al. [2009]. The introduction of gas in our model would likely increase surface velocities by several orders of magnitude compared with the present simulations, which are controlled solely by crystal settling and thermal convection. A higher frequency of pulsation of these surface velocities can also be expected because gas motion will affect the distribution of surface temperatures, which will result in non-uniform surface velocities. The release of gas is also likely to increase the rate of crystallization; however the latent heat of crystallization would act as a buffer that would compensate the cooling induced by this release of gas to an extent that we have not quantified.

[76] We assumed that the Erebus system contained already formed crystals that are distributed homogeneously and in equilibrium with the melt at all times. Considering that the simulations started from an assumed constant crystallinity (30 vol.%), our simulations show that between 6 and 20 vol.% crystals remain in suspension, depending on whether the system is closed or open, respectively. Crystal content observed in phonolite lavas from summit (∼30–40 vol.%) and lava bombs (∼30 vol.%) [Kelly et al., 2008] are the product of a thermal and pressure history inside the magma system. If those natural products are indicative of equilibrium conditions within the system, then our calculated values are below the observed amount of crystals and the differences are larger in closed system than in open system. This can be explained by the fact that closed system simulations compensate cooling-induced contraction by underestimating crystal content. A more quantitative assessment is difficult; one reason is that the extrapolation toward low pressures and low water content of existent phase equilibrium data on phonolitic magmas [Andújar et al., 2008, 2010] is subject to caution, notably because the fluid composition (H₂O, F, Cl and Fe) may have a strong effect on the presence or absence of phases and their compositions. Crystallization may play an important role to fill the gap between our results and natural observations; however the theory developed by Lavorel and Le Bars [2009] and our simulation equally indicate that it will compete with particle settling and convective transport. Crystal sizes in Erebus lava suggest that crystallization at shallow levels occurs more by growth of pre-existing crystals than by addition of newly nucleated crystals. We thus anticipate that improved knowledge of low-pressure phonolite petrology and feldspar growth kinetics could bring crystallization laws that, once included in our simulations, may resolve these differences in crystal volume fraction. We nevertheless posit that the main observations and conclusions developed in this paper would remain largely unchanged for steady state simulations including crystallization.

[77] The unusually large crystal size at Erebus controls not only when convection changes from being thermally driven to being settling driven, but also how much of the crystal load remains in suspension. Crystals in other lava lakes, such as Nyiragongo [Demant et al., 1994], Villarrica [Witter et al., 2004] or Pu' u O' [Garcia et al., 2000] are one to two orders of magnitude smaller in size. Their influence on the convective pattern is thus smaller than at Erebus. Assuming crystal sizes of 5 mm and using equation (13), Erebus lava lake would have an initial value but a value at equilibrium N_eq = 66000 particles/cell with an efficiency factor of one. Lowering N_eq to reach a solution where 100% of crystals would remain in suspension suppose a very low value of the efficiency factor (σ = 10⁻²), far away from the order of 1 usually accepted for this constant. We instead simply suggest qualitatively that means that energy balance would be in favor of convective motions and that most 5 mm crystals would remain in suspension. We can thus infer that the convection within lava lakes bearing crystals ≪1 cm is quickly thermally dominated, and that crystal size is an important control of how crystal-rich the magma may remain, which, in turn, controls suspension rheology.

6. Conclusions

[78] We have studied the convective patterns beneath the surface of a lava lake by means of numerical simulations. The bi-phase, fluid dynamical model was validated by comparing the simulated structures of convective currents with those measured in the analogue experiment designed to represent cooling magma performed by Jaupart and Brandeis [1986]. Cooling generates Rayleigh-Taylor instabilities, the simulated motions of which are consistent to first order with those produced experimentally. The velocity at which the modeled plumes sink are in good agreement with that predicted by theoretical models and differs only within the bounds of experimental uncertainty from that seen in the experiment. At longer times, the temperature evolution in the well-mixed, central part of the tank is consistent with that observed in the experiments, though time shifted because of inevitable differences in initial conditions. The bi-phase capability of the numerical model was tested by calculating the bulk viscosity of a series of Poiseuille flow simulations with various amounts of particles. Results suggest that the model is applicable to a mixture of magmatic liquid and crystals as long as the crystal load is at or below close packing (65–67 vol.% crystals in this work).

[79] We applied the validated and verified model to an idealized Erebus magmatic system composed of a lava lake, a feeding conduit, and a magma chamber. Model parameterization was as close to the natural system as allowed by model assumptions. We tested the influence of conduit geometry, crystal load, and the nature of system boundaries on the convective regime. Overall, model outputs show that the steady state convective activity is sufficient to keep the magma above glass transition throughout a 30 yr period provided that the conduit diameter is large enough (10-m diameter), but regardless of whether the system is filled with a pure melt, a perfectly homogeneous mixture of magma with 30 vol.% crystals, or a bi-phase magma with an initial content of 30 vol.% crystals. If the conduit is narrow (4 m), the central viscous instability descending within the lake from the quenched surface chokes off the conduit and prevents hot fluid from ascending from the chamber, which slows down the whole convection rate. This phenomenon occurs at different times for the mixture and bi-phase simulations, but never for the pure melt because these simulations fail to reach steady state before the considered 30 years.

[80] Combining together crystals and melt into a virtual single-phase can reproduce some features of such systems, such as (i) the early formation of a central instability in the lake, (ii) the global temperature evolution and (iii) the order of magnitude of the average velocity range of the flow motion at the surface of the lake. The strength of the bi-phase approach, on the other hand, is to take crystal settling into account. For instance, a series of changing convective patterns are observed for bi-phase simulations throughout the simulated period. We found that conduit diameter strongly influences the circulation within the conduit itself. Overall, conduit- and system-wide patterns are in good agreement with different experiments and related numerical models. Our results show that all these emerging patterns can be gathered into a single framework. This very richness of behavior, however, poses a challenge to design numerical simulations with testable model outputs.

[81] Crystals are efficiently transported by the liquid phase but a small decoupling due to their large size (5 cm) causes settling. Compared to a crystal-free magma, their presence accelerates the convective rate and changes the convective regime, enhancing the effective heat transfer between lake and chamber. Thermal forces initially prevail and the suspended crystal content decreases following a power law until the particles accumulate in the magma chamber without being re-circulated in the whole system. A settling-driven convection then dominates, during which between 6 and 20 vol.% crystals remain in suspension, depending on whether the system is closed or open, respectively. Compared to closed-system simulations, an open system results in higher temperatures in the whole domain, higher convective rates (in both frequency and amplitude) at the surface of the lake, and a higher suspended crystal fraction, closer to the initial crystallinity. We observe a lasting crystal-rich layer of 10 m at the bottom of the chamber. Closed system simulations compensate cooling-induced contraction by underestimating crystal content. As a result, open system simulations should be more realistic than closed system ones because they correctly account for cooling-induced volume reduction.

[82] Our stringent assumptions lead to an idealized system with rich dynamics but that only poorly compares with natural observations. The average velocity at the lake surface is on the order of 10⁻⁶ m/s (closed system) and 10⁻⁵ m/s (open system). These simulated crystal loads and lake surface velocities in steady state are much lower than those observed at Erebus (30 vol.% and 10⁻¹ m/s, respectively). The permanent release of gas at Erebus, unaccounted for in our simulations, is thus likely to contribute greatly in sustaining such a vigorous convection, preventing crystal settling and increasing the velocity of magma at the surface. We nevertheless estimate that the unusually large crystal size seen in Erebus phonolite controls not only when the convection changes from being thermally driven to being settling driven, but also how much of the crystal load remains in suspension. This is probably not true for other magmatic systems with much smaller crystals.