Articles | Volume 15, issue 9
Geosci. Model Dev., 15, 3773–3796, 2022
Geosci. Model Dev., 15, 3773–3796, 2022
Model description paper
10 May 2022
Model description paper | 10 May 2022

MagmaFOAM-1.0: a modular framework for the simulation of magmatic systems

MagmaFOAM-1.0: a modular framework for the simulation of magmatic systems
Federico Brogi1,2, Simone Colucci1, Jacopo Matrone3, Chiara Paola Montagna1, Mattia De' Michieli Vitturi1,4, and Paolo Papale1 Federico Brogi et al.
  • 1Istituto Nazionale di Geofisica e Vulcanologia, sezione di Pisa, Pisa, Italy
  • 2Istituto Nazionale di Oceanografia e di Geofisica Sperimentale, Trieste, Italy
  • 3Dipartimento di Matematica, Università degli Studi di Firenze, Florence, Italy
  • 4Department of Geology, University at Buffalo, Buffalo, New York, USA

Correspondence: Federico Brogi (


Numerical simulations of volcanic processes play a fundamental role in understanding the dynamics of magma storage, ascent, and eruption. The recent extraordinary progress in computer performance and improvements in numerical modeling techniques allow simulating multiphase systems in mechanical and thermodynamical disequilibrium. Nonetheless, the growing complexity of these simulations requires the development of flexible computational tools that can easily switch between sub-models and solution techniques. In this work we present MagmaFOAM, a library based on the open-source computational fluid dynamics software OpenFOAM that incorporates models for solving the dynamics of multiphase, multicomponent magmatic systems. Retaining the modular structure of OpenFOAM, MagmaFOAM allows runtime selection of the solution technique depending on the physics of the specific process and sets a solid framework for in-house and community model development, testing, and comparison. MagmaFOAM models thermomechanical nonequilibrium phase coupling and phase change, and it implements state-of-the-art multiple volatile saturation models and constitutive equations with composition-dependent and space–time local computation of thermodynamic and transport properties. Code testing is performed using different multiphase modeling approaches for processes relevant to magmatic systems: Rayleigh–Taylor instability for buoyancy-driven magmatic processes, multiphase shock tube simulations propaedeutical to conduit dynamics studies, and bubble growth and breakage in basaltic melts. Benchmark simulations illustrate the capabilities and potential of MagmaFOAM to account for the variety of nonlinear physical and thermodynamical processes characterizing the dynamics of volcanic systems.

1 Introduction

Simulating transport processes in volcanic systems is of crucial importance to understand the physics of eruptions, correctly interpret geophysical signals recorded by volcano monitoring systems, anticipate volcanic scenarios, and forecast volcanic hazards (Sparks2003; Bagagli et al.2017). A great number of flow models have been developed to address specific volcanic processes, including magma chamber dynamics (Ruprecht et al.2008; Bergantz et al.2015; Garg et al.2019), conduit flow (Melnik2000; Papale2001; de' Michieli Vitturi et al.2008b; Colucci et al.2017b), volcanic plumes (Suzuki et al.2005; Cerminara et al.2016), pyroclastic flows (Esposti Ongaro et al.2007; de' Michieli Vitturi et al.2015; Dufek2016), and lava flows (Griffiths2000). Inter-model comparison studies have evaluated individual model performance and the relevance of the different subprocesses, and they have highlighted target areas for improvement (Massol and Koyaguchi2005; Macedonio et al.2005; Sahagian and Proussevitch2005; Costa et al.2016). All these models attempt to tackle the great complexity arising from the presence of multiple phases. Interactions among liquid phases (e.g., silicate melt), solid phases (e.g., crystals or pyroclasts), and gas phases (exsolved volatiles or atmospheric gas) are indeed ubiquitous in volcanic systems, from deep magma chambers up into the atmosphere (e.g., Jackson et al.2018; Keller and Suckale2019).

Volcanic transport processes are typically characterized by a wide range of spatial and temporal scales at which different interacting physical subprocesses occur (Griffiths2000; Gonnermann and Manga2007; Dufek2016). From a modeling perspective, there is no general approach able to treat all these subprocesses at the same time, and thus specific models are usually developed for each application.

Figure 1Schematic representation of the (a) multiscale nature of a gas–liquid–solid flow (reprinted with minor modifications from Yeoh and Tu2019b, with permission from Elsevier), (b) a two-phase volume at the local scale containing discrete phase constituents (e.g., bubbles or droplets), and the corresponding approximation path to describe it as made of two interpenetrating continuum fluids (modified from Marschall2011).

A generic multiphase system can be thought of as composed of sub-domains or regions pertaining to single phases (gas, liquid, or solid) separated by interfaces (boundaries) representing sharp discontinuities where the physical properties change abruptly. The main challenge in modeling multiphase with respect to single-phase flows is due to the presence of such a discontinuity. The topology of this interface defines the amount of interfacial area that is available for the phases to exchange mass, momentum, and energy and strongly affects the behavior of the multiphase mixture. Moreover, this interface is not static but changes dynamically with the flow, and complex flow features may emerge due to the presence of moving phase boundaries. Understanding and modeling multiphase flows also requires taking appropriate consideration of their multiscale character. The typical size of the interfaces can be comparable to, or orders of magnitude smaller than, the domain and flow length scales, or it can even cover a broad range of scales (Fig. 1). Cascading effects and multiple coupled phenomena at the different scales may have dramatic consequences for the flow dynamics. At the scale of the system (e.g., volcanic conduit), large flow structures that govern the flow directly depend on the properties of the multiphase mixture, which in turn are determined by the dynamic reorganization of local mesoscale structures (e.g., coalescence and/or breakage of large bubbles or bubble cluster dynamics) as well as the motion of individual constituents at the microscale (e.g., single small bubbles) within a continuum phase (e.g., liquid). Modeling such complex phase interactions across a wide range of scales certainly represents a big challenge.

Interface-resolving methods, similar to direct numerical simulation (DNS) approaches in single-phase turbulent flows (Moin and Mahesh1998), fully resolve the scales of the fluid equations and track the topology of the interfaces. With this approach no assumptions are made regarding the properties of the multiphase mixture or interfacial phase exchanges. The dynamics of the multiphase system emerge naturally from the computation as a direct consequence of solving phase interactions locally at the interfacial scale (under the constraints of mass, momentum, and energy conservation for each phase). DNS can therefore be thought of as a virtual laboratory to understand fundamental physics, especially at the microscale, that can capture emerging dynamics resulting from nonlinear phase interactions that are difficult to parameterize a priori (e.g., Segre et al.2001). DNS can also provide a detailed description of the flow that often is not accessible in experiments. Therefore, it can help in interpreting laboratory observations (e.g., Qin and Suckale2020) or even building and testing constitutive–parametric models for interphase interaction exchanges and the behavior of the multiphase mixture as a whole (Fang et al.2019). In volcanology, the DNS approach has been used to study large gas bubbles ascending in a conduit through low-viscosity melts (Suckale et al.2010a), buoyancy-driven instabilities in liquids at different densities (Suckale et al.2010b), and the complex rheological behavior of crystal-bearing magmas (Qin and Suckale2020). Based on the computationally efficient lattice Boltzmann method, interface-resolving modeling has also been useful to better understand bubble growth, deformation, and coalescence (Huber et al.2014) as well as the mush microphysics characterizing crystal-rich magma reservoirs (Parmigiani et al.2014). Despite its proven ability to provide understanding of the fundamental physics near interfaces, DNS remains limited to specific, computationally tractable problems, since it requires large of computational resources.

Simulations of magmatic systems that can aid the interpretation of geophysical (Bagagli et al.2017) or petrological observations (Cheng et al.2020) need to cope with very large domains on the scale of the kilometers. On the other hand, the smallest scales (e.g., small eddies, bubbles, crystals, Fig. 1a) remain several orders of magnitudes smaller, resulting in skyrocketing computational costs for DNS even when considering relatively simple flow conditions (e.g., laminar flows, Yeoh and Tu2019a). Multiphase flows that present dispersed interfaces (e.g., bubbly, droplet or particle-laden flows, Fig. 1b), which are common in both natural and industrial settings (Keller and Suckale2019; Moreno Soto et al.2019), have been successfully modeled with a different approach. The multi-fluid formulation employs averaging techniques that filter out the interfacial scales that are too small to be resolved (Marschall2011). The complexity of a volume with multiple phases at the local scale is characterized by phase-average properties and a volumetric fraction that expresses the relative presence of one phase with respect to the others (Fig. 1b). Neglecting the details of the topology of the interfaces at the local scale allows describing the phases at the system scale as interpenetrating continua governed by separate sets of conservation equations. The resulting equations hence resemble those for single-phase flows except for the volumetric fraction and the presence of phase interaction terms that require appropriate closure. Similarly to large eddy simulations for turbulent flows, additional constitutive models are in fact required to recover the physics of the missing small scales. The multi-fluid approach allows modeling thermomechanical disequilibrium (e.g., phases with different velocities, temperatures, or pressures) as well as interactions of the dispersed phases for any multiphase system (Marchisio and Fox2007), including magmas (Keller and Suckale2019). Applications of multi-fluid modeling in volcanology include but are not limited to the study of buoyancy-driven magma mixing (Ruprecht et al.2008), conduit dynamics (Papale2001; Dufek and Bergantz2005), and volcanic plumes (Neri et al.2003; Ongaro et al.2007). However, the definition of appropriate closure models for interfacial phase exchanges is crucial for these models to provide accurate predictions of the evolution of the multiphase system. The closure of the multi-fluid formulation remains challenging and is mostly achieved using system-dependent constitutive equations that are often empirical and valid for specific flow regimes (in terms of interface topology and/or concentrations of the dispersed phase). Thus, the generality of the multi-fluid formulation is reduced by the specificity of the constitutive models. Multi-fluid models are, however, more computationally expensive than single-phase models, as they require an additional set of governing equations for each phase. As the number of phases increases, the computational burden also increases dramatically (e.g., Ferry and Balachandar2001). The definition of the interfacial exchange terms can also indirectly increase the computational cost. For instance, the fluid–particle drag introduces a timescale in the equations, which is the relaxation time of the dispersed phase, that describes the time required by the particle to adapt to a change in velocity of the surrounding fluid. When this relaxation time is small (typically for small particles and/or high fluid viscosities), the stability and accuracy of the numerical solution require a time step smaller than the relaxation time, increasing the number of iterations needed to solve the flow timescale. Under the assumption of thermomechanical equilibrium the equations of the multi-fluid model can be further reduced to an evolution equation for a single pseudo-fluid representing a mixture of multiple phases. From a computational point of view, given the reduced number of equations needed to track the evolution of the mixture, this is a more convenient approach. In addition, when there is a strong thermomechanical coupling between phases (small relaxation times), it is reasonable to assume that the particle velocity is equal to the fluid velocity, effectively removing the aforementioned issues related to the definition of the interaction terms and the relaxation time. The single mixture, pseudo-fluid approach has received some success within the volcanological community. Simulations of magma mixing (Longo et al.2012; Garg et al.2019) and conduit dynamics (de' Michieli Vitturi et al.2008a; Melnik and Sparks2006) as well as volcanic plumes (e.g., Suzuki et al.2005) are only a few examples of application. Nevertheless, as for the multi-fluid, constitutive models play a crucial role in determining the reliability of single-fluid predictions. Finally, it is also important to emphasize that both the multi-fluid and the single-fluid mixture models are based on the average form of the multiphase flow equations. Averaging procedures implicitly require the separation of scales. In volume averaging, for instance, the volume should be large enough to contain a representative sample of the dispersed phase (e.g., bubbles or particles) but much smaller than the typical distance over which flow properties vary significantly. This condition is rarely satisfied in real applications since intermediate scales between the local and the system scale are present (Brennen and Brennen2005). In spite of all the issues and limitations, a few of which have been mentioned above, a large amount of theoretical and practical investigations remain based on single- and multi-fluid models (Yeoh and Tu2019a).

The increased ability of models to include detailed physics strictly requires the development of more flexible computational tools that can easily switch between constitutive models and solution techniques to adapt to different dynamical regimes, thereby reducing computational efforts, increasing usability, and easily allowing scientists to perform inter-model comparison studies and models coupling.

The open-source library OpenFOAM provides a variety of fluid solvers for multiphase flows that can be combined with several different constitutive equations. Its modular object-oriented implementation allows developers to easily expand and adapt the code and users to combine different models at runtime with almost no need to code. Given a set of discretized fluid evolution equations (or “solver”), the user can easily select appropriate thermophysical and rheological models or switch from 2D or 3D to axis or plane symmetric simulations. The OpenFOAM community is continuously growing, as is the range of applications of interest for both the academy and industry (e.g., Winden2021). Moreover, the recently established exaFoam consortium will improve computational performance, enabling the “OpenFOAM community” to efficiently exploit “the current evolving HPC hardware and middleware” (, last access: 27 April 2022). OpenFOAM is thus an ideal platform for developing a computational toolbox for the next generation of magmatic systems modeling. In this work we present the MagmaFOAM library, an extension of OpenFOAM dedicated to solving multiphase volcanic flows. The current implementation features multiple volatile saturation models (Papale et al.2006) and specific formulations for the equation of state (Lange and Carmichael1987) and viscosity (Giordano et al.2008) of magmatic mixtures including dissolved volatiles. MagmaFOAM retains the basic coding principles of OpenFOAM, inherits its flexibility, and takes full advantage of the family of fluid solvers and constitutive models (e.g., non-Newtonian rheological models) already implemented in OpenFOAM.

This paper is structured as follows. First we provide an overview of the basic ingredients of MagmaFOAM, including the specific magmatic constitutive equations and how they are implemented. Then, we show benchmarks and validation tests aimed at verifying the code ability to solve problems for segregated and dispersed flows of interest for magmatic systems with different modeling approaches. Finally, we summarize and discuss our results and draw the conclusions.

2 MagmaFOAM ingredients

2.1 Structure of MagmaFOAM

MagmaFOAM uses the same organization of OpenFOAM (Fig. 2), and its hierarchy is therefore subdivided into applications and libraries (src). Code organization is therefore rational and efficient, reducing code duplication, promoting code re-usage, and facilitating testing. Most of the applications are assembled at runtime based on user requests using dynamic linking to pre-compiled libraries: before running a simulation the user can arbitrarily select boundary conditions, discretization schemes, and mixture and phase constitutive equations. This mechanism allows selecting and combining modeling ingredients, among the possible combinations, from both OpenFOAM and MagmaFOAM (Fig. 2) without the need for coding.

Figure 2MagmaFOAM–OpenFOAM coupling scheme.


Multicomponent constitutive models for magmatic systems

The dynamics of magmas as they ascend, stall through the crust, and possibly erupt are strongly dependent on their physical properties (mostly density ρ and viscosity μ), which in turn are determined by composition and phase distribution as well as pressure p and temperature T conditions. The interplay among pT conditions, melt–crystals–bubbles phase changes, and density and viscosity variations originates a wealth of possible space and time patterns for magma storage, transport, and eruption (e.g., Lesher and Spera2015). When handling this thermophysical complexity, state-of-the-art multicomponent constitutive models that compute melt properties as a function of the local pressure, temperature, and composition are the necessary basic ingredients and have been implemented in MagmaFOAM.

Multicomponent volatile saturation is included through the SOLWCAD model (Papale et al.2006), which provides equilibrium H2OCO2 saturation over a broad range of pT conditions and for virtually any melt composition. This model overcomes the ideal Henrian behavior, which is a reasonable approximation only at low pressures (100 MPa; Papale et al.2006). Once the phase distribution of volatile species is computed through SOLWCAD, the relevant physical properties for the multiphase magma can be derived.

The density of the silicate melt up to a few gigapascals (GPa) (3 GPa or 100 km depth) is computed as in Lange and Carmichael (1987) with an empirical equation of state as a ratio of the oxides' molar masses (Mi) and molar volumes (Vi):

(1) ρ ( p , T , X ) = M ( X ) V ( p , T , X ) = i X i M i i X i V i ( p , T ) ,

where Xi is the mole fraction of the ith oxide component. To a good approximation, molar volumes do not depend on melt composition (Lesher and Spera2015) and can be computed with a polynomial expansion:

(2) V i ( p , T ) = l , m a l , m i T l p m = a 0 , 0 i + a 1 , 0 i T + a 0 , 1 i p + a 1 , 1 i p T +

The polynomial coefficients al,mi have been determined from laboratory experiments. For the oxides we have used the coefficients reported by Lange and Carmichael (1987) and Lesher and Spera (2015). For H2O and CO2 we referred to Burnham and Davis (1974) and Papale (1999), respectively.

Melt viscosity is described as in Giordano et al. (2008). This model includes temperature and compositional effects for a wide range of melt compositions. In addition, the model can be used to determine the compositional dependence of important viscosity-derived properties, such as melt glass transition temperature and fragility. This aspect may be particularly relevant when modeling the ascent of degassing magma to determine the potential for brittle fragmentation. A drawback is that the model does not take into account the effect of pressure on viscosity, which can become relevant when modeling magma transport in the deep crust and mantle.

The model is based on the Tammann–Vogel–Fulcher (TFV) relationship for the non-Arrhenian temperature dependence of the bubble–crystal free viscosity η:

(3) log η = A + B T - C ,

where T [K] is the temperature, A is a constant, and B and C are parameters that depend on the melt composition, including dissolved volatile species. The A constant provides a high-temperature limit for viscosity (10-4 Pa s) that holds for all melts regardless of their composition and is supported by both theoretical considerations and experimental observations (e.g., Scopigno et al.2003). Let us also note that Eq. (3) has no physical meaning for TC (Mauro et al.2009).

2.2 Modeling volatiles concentration at the bubble–melt interface

Models accounting for multicomponent phase change require a description of the evolution of the composition at the interface between phases. The mass transfer rate (per unit volume of liquid + gas) of a volatile component can be defined as the product between the interfacial mass flux Ji [kg (m2 s)−1] and the interfacial area concentration A [m2 m−3].

(4) Γ i = J i A

The area concentration A is determined by the geometrical configuration of the gas–liquid interface, and hence it is strongly dependent on the flow regime. It can be computed using simple geometrical assumptions on the dispersed phase (e.g., monodisperse bubbles with constant radius) or, for more complex flow scenarios, with additional transport equations (e.g., IATE model; Ishi and Hibiki2006). The model for Ji expresses the driving force for diffusive mass transfer of the component i and can be calculated with the following relationship:

(5) J i = k i Δ C i ,

where ki [kg (m2 s)−1] is the mass transfer coefficient, a function of the diffusion coefficient Di (Cussler2009; Thummala2016), and ΔCi is the difference between the mass fraction of the species in the phase (Ci) and at the interface (Cif):

(6) Δ C i = C i - C i f .

Under the assumption of local equilibrium, the mass fraction at the interface can be expressed as

(7) C i f = C i sat p , T f , X , X V tot ,

where Cisat is the saturation concentration of a specific volatile species (i.e., mass fraction at thermodynamic equilibrium). In general this is a nonlinear function of pressure (p), temperature at the interface,(Tf), melt composition (X), and total amount of volatiles of all species (XVtot). For magmas with H2O and CO2, Cisat can be computed using SOLWCAD (Papale et al.2006) or other dedicated models (e.g., Newman and Lowenstern2002; Burgisser et al.2015). Direct coupling of any fluid solver with these models is usually too computationally expensive. Therefore, MagmaFOAM solvers can read the saturation surface from a pre-processed table. During the simulation, tabulate values are interpolated (multilinear interpolation) and used to compute Cisat in Eq. (7).

2.3 MagmaFOAM constitutive models

Constitutive models implemented in MagmaFOAM can be selected and combined at runtime (no need for coding) with existing OpenFOAM solvers suitable for the specific problem under consideration (Fig. 2). For example, the MagmaFOAM model for silicate melt density can be used with any compressible solver, either single-phase or multiphase. This constitutive model is not compatible with incompressible solvers that require density to be constant; however, in this case the density of the incompressible fluid can be preliminarily defined by taking advantage of the dedicated MagmaFOAM utility magmaThermoMixture. The latter can also be used for testing implemented models as it simply returns the thermophysical properties as a function of composition, pressure, and temperature. Demonstrative tutorials are included in MagmaFOAM to show how the end user can accomplish all these tasks at runtime using both single-phase and multiphase solvers.

2.4 Models for multicomponent bubble growth

Volatile phase changes and bubble growth are ubiquitous processes in volcano dynamics (Proussevitch and Sahagian1998). The gas exsolution process begins with the nucleation of bubbles in an oversaturated melt and continues with bubble growth. Bubbles grow by mass diffusion when the silicate melt is oversaturated in volatiles and by mechanical expansion as a response to pressure decrease. The viscosity of the surrounding melt and the surface tension oppose a resistance to bubble growth and control the mechanical disequilibrium between the bubbles and the melt itself. A number of works (Proussevitch et al.1993; Lyakhovsky et al.1996; Proussevitch and Sahagian1998; Lensky et al.2001, 2004; Chouet et al.2006; Shimomura et al.2006; Coumans et al.2020) solve the system of bubbles as a monodisperse periodic array of static, spherical, single-component (H2O) growing bubbles surrounded by a viscous melt shell using the Rayleigh–Plesset equation. A suite of models based on a similar approach have been implemented in MagmaFOAM and benchmarked to simulate multicomponent diffusive bubble growth. This approach provides, at low computational cost, an accurate representation of the coupled momentum balance and diffusive transport of volatiles because it resolves the concentration profile well near the bubble interface (Huber et al.2014). The strong assumptions that the size distribution is monodisperse and the bubbles are non-deformable and mechanically coupled with melt limit the range of applicability of the model. In high-viscosity systems at low vesicularity, the model can provide reliable results when compared with experiments (e.g., Coumans et al.2020). The model does not take into account interfacial interactions (fluid–particle and particle–particle) that can give rise to emergent behavior, as in the case, for example, of bubble waves (Manga1996). All model equations can be found in Appendix B and are solved as systems of ordinary differential equations (ODEs) using the OpenFOAM ODE solvers.

3 Benchmarks and test cases

The test cases presented here are included in the MagmaFOAM distribution together with the relevant post-processing routines. The results shown here are thus fully reproducible, and the benchmarks can be used to study the accuracy and efficiency of other OpenFOAM or external solvers.

3.1 Interface-resolving modeling

The volume of fluid (VOF) method is adopted in OpenFOAM to resolve the position and shape of the interface separating two fluids or phases (e.g., liquid–gas). This methodology treats the interface discontinuity as a smooth but rapid variation (few computational cells) of an indicator field (volumetric fraction) representing the relative presence of one phase with respect to the other in each cell. The volumetric fraction is zero or 1 away from the interface, allowing distinguishing between one phase and the other, and assumes intermediate values in the region containing the interface. As a result, the location of the interface and its shape are known only implicitly from the volumetric fraction. The evolution of the interface is then obtained by simply advecting the volumetric fraction using the velocity field computed from a single-fluid (e.g., the OpenFOAM solver interFoam) or multi-fluid momentum equation (e.g., the OpenFOAM solver multiphaseEulerFoam). The transport equation for the indicator function is under the constraint of mass conservation, and therefore, with respect to other methods (e.g., level-set method), VOF is mass-conservative by construction. However, in practice the conservation of mass depends on the accuracy in numerically solving the transport equation. The discontinuous nature of the volumetric fraction (a step function) at the interface makes the numerical solution of this equation challenging. In particular, numerical diffusion due to the discretization of the advection term prevents a sharp representation of the interface that tends to be smeared over the computational cells, causing inaccurate estimations of its position and curvature. Different techniques exist to solve this issue. With a geometrical approach one may reconstruct the position of the discontinuity at the subgrid level, provided that the interface can be described with a specific functional form (Rider and Kothe1998; Aulisa et al.2003). The interface is then advected by the flow in a Lagrangian manner. This technique effectively prevents numerical diffusion and provides a more accurate representation of the interface at the cost of a significantly more complex algorithm and increased computational load. Other approaches rely on relatively more simple algebraic solutions that reduce numerical diffusion (e.g., Ubbink and Issa1999). Specifically, interFoam makes use of a high-order differencing scheme (in the interface region only) and an additional compressive term in the advection equation that effectively counterbalances the numerical diffusion of the interface. While this approach is simpler and less computationally expensive than the geometrical reconstruction, the interface is spread over few computational cells and its precise position remains unknown. Nevertheless, in kinematic tests, interFoam has shown good mass conservation properties and acceptable advection errors (Deshpande et al.2012). Spurious currents and artificial deformations of the interface are also an issue with VOF. Inaccurate interface curvature, together with a discrete force imbalance at the interface, typically produces spurious vortices that can artificially deform the interface. Depending on the simulation setup, the kinetic energy of these vortices may rapidly decay or grow and in the worst case scenario even cause the simulation to crash. However, spurious currents may pose a serious issue, mostly for surface-tension-dominated flows, and are less important for inertia-dominated flows. For interFoam, Deshpande et al. (2012) have shown that the growth of spurious currents can be controlled by choosing an appropriate time step; interFoam solves flows characterized by constant, or slowly varying with respect to the flow timescales, fluid properties. This approximation holds for relevant volcanic scenarios as gas-poor magmatic reservoirs at depth, which are characterized by relatively fast overturn times (Ruprecht et al.2008; Perugini et al.2010; Montagna et al.2015).

Here we present benchmarks and test cases to evaluate the accuracy of the solver interFoam to explore the dynamics of two immiscible fluids separated by a free interface. Specifically, we perform detailed studies of buoyancy-driven magma mixing and rising bubble dynamics. Overall, we find remarkably good agreement between our simulation results and theoretical or numerical results from the literature over different flow regimes of interest for magma dynamics. The numerical solutions relative to cases with low Reynolds numbers Re are very accurate (e.g., Figs. 4 and 7). At larger Re, the results are less accurate due to the appearance of high-frequency numerical noise that can trigger secondary spurious interface instabilities (e.g., Fig. 5). Reducing numerical noise by adopting different numerical schemes is one relevant element for future investigation. The magnitude of the compressive term, used in the solver to prevent numerical smearing of the interface, is a free parameter in the simulations and may influence the accuracy of the solution depending on the problem parameters. More recent OpenFOAM versions include more rigorous and accurate interface-resolving methods (e.g., Roenby et al.2017).

3.1.1 Interacting magmas

Magma is thought to rise from the mantle into the crust in discrete batches (Annen et al.2006) that then tend to stall and cool at different depths, while their chemistry evolves towards more felsic compositions (Sigurdsson et al.2015). Different batches of magma may interact as they ascend towards shallower depths, resulting in magma mingling and mixing. The latter are widespread phenomena in volcanic plumbing systems (Perugini and Poli2012; Morgavi et al.2017) and have often been invoked as eruption triggers (Wark et al.2007; Druitt et al.2012; Martí et al.2020). Mingling and mixing are typically driven either by gravitational Rayleigh–Taylor instabilities, involving contacts between magmas with different densities due to compositional, thermal, or phase stratifications (e.g., Jellinek et al.1999; Montagna et al.2015; Garg et al.2019), or by percolation of pressurized magmas arriving from depth into mushy reservoirs (Bachmann and Bergantz2003; Seropian et al.2018).

A standard benchmark to test numerical solvers for Rayleigh–Taylor instability problems requires comparing computed growth rates for small-amplitude single-mode perturbations with the linear stability theory. The latter predicts that a small perturbation grows exponentially with a rate that depends on its wavelength, fluid density and viscosity contrasts (Chandrasekhar1955), surface tension (Chandrasekhar2013), compressibility (Mitchner and Landshoff1964), and diffusivity (Duff et al.1962; Xie et al.2017). The problem parameters can be expressed by two dimensionless numbers: the Atwood number Atw=(ρh-ρl)/(ρh+ρl) and the Reynolds number (Re=WgW/ν), where ρh and ρl are the two liquid densities, ν is the kinematic viscosity (νh=νl), W is the wavelength of the perturbation, and g is the gravitational acceleration. We consider a 2D rectangular domain with a no-slip condition (walls) on top and bottom boundaries and periodic conditions on the sides. The interface between the two liquids is located at the center of the computational domain (Fig. 4). The size L of the computational box is determined by the wavelength of the initial perturbation (L=W×2W). Benchmark results are reported in Fig. 3 for Atwood numbers relevant for natural melts. The computed growth rates are in agreement with the theory (Xie et al.2017) for different wavelengths (or equivalently wave numbers k=2π/W) of the perturbation. The solver underestimates the peak growth rates at low k, corresponding to high Re. A more in-depth analysis of the results (Appendix A) reveals that this discrepancy is mainly due to an initial delay in the onset of the perturbation. Removing this initial offset, the computed growth rates are much more accurate. Smaller initial perturbation amplitudes also improve accuracy.

Figure 3Comparison between computed growth rates (symbols) of the Rayleigh–Taylor instability in the linear regime obtained with the solver interFoam and theoretical ones (dashed line). Bubble growth rates are computed by tracking the position of the interface with respect to the central axis of the domain, while spike growth rates are computed with respect to one of the lateral boundaries.


Figure 4Rayleigh–Taylor instability (A=0.5,Re=256) computed with the OpenFOAM solver interFoam. The density field (color-coded) is compared with the density contours in He et al. (1999) (black lines).

Figure 5Rayleigh–Taylor instability (A=0.5, Re=2048) computed with the OpenFOAM solver interFoam (in color) with (a) 256×1024 cells and the default value for the interface compression factor Cα=1 as well as (b) 512×2048 cells and Cα=0.1. The density field (color-coded) is compared with the density contours of He et al. (1999) (black lines).

As the instability grows and its amplitude becomes comparable with its wavelength, nonlinear effects become dominant and the linear theory is no longer valid to predict the evolution of the system. In order to validate interFoam for nonlinear regimes we have compared our results with He et al. (1999) for single-mode perturbation with a 10 % amplitude-to-wavelength ratio, Atwood number A=0.5, and Reynolds number Re=256. Remarkably good agreement is obtained for the evolution of the fluid interface (Fig. 4) using the same resolution (256×1024 cells). Convergence of the results was tested for different space and time resolutions using the adaptive time step based on the maximum allowed Courant number (Comax=0.5) to speed up the simulation. For a given mesh resolution the accuracy and convergence of the solution depend on the values of Comax and the number of iterations (nIter) used to solve the pressure–velocity coupling with the PISO algorithm (Issa1986). Generally, larger Comax (<1 for numerical stability) values require larger nIter for solution convergence; our experience suggests that a relatively high number of nIter balances larger values for Comax, reducing computational times. This way, even if the errors relative to the continuity equation are larger, the solution is not significantly affected.

For the high-Reynolds-number test case (Re=2048) of He et al. (1999), the quality of the solution deteriorates using the same resolution (Fig. 5). The interface is deformed by artificial secondary instabilities most probably triggered by spurious numerical noise. Removing the interface compression term and doubling the number of cells improves the solution to nearly match the reference.

Magmas usually interact both mechanically and chemically, and therefore the immiscible approximation described above is not justified a priori. Nevertheless, to first approximation and on relatively short timescales (hours to days), chemical diffusion among interacting magmas at the plumbing system scale can be neglected (e.g., Ruprecht et al.2008; Garg et al.2019), and magmas can be considered immiscible. Here we describe exemplary buoyancy-driven interaction among two natural silicate melts (Fig. 6). Density and viscosity of the two melts are computed a priori using the MagmaFOAM utility Test-magmaThermoMixture. As a test case, we reproduce at small scale a typical (Garg et al.2019) interaction among a volatile-rich basalt (XH2O=2 wt %) and a chemically more evolved andesitic melt. Temperature is set to T=1300C and pressure is atmospheric. Melt compositions are reported in Table D1. The composition as well as p and T conditions are considered only in the pre-processing to compute the density and viscosity of the melt that remain constant throughout the simulation. The relevant dimensionless numbers are now A=0.0167982 and Re=54.065 for a physical domain 1 m×4 m. Surface tension is again neglected. Compared to the previous simulations (e.g., Fig. 4), the two liquids now have different kinematic viscosities. The larger viscosity ratio requires significantly increasing the numbers of iterations (≈300) needed to solve pressure–velocity coupling (keeping Comax=0.5). As a result, the simulation is computationally much more demanding. The simulated time covers the entire overturning process (Fig. 6).

Figure 6Rayleigh–Taylor instability between a volatile-rich (XH2O=2 wt %) basalt (bottom) and andesite (top) computed with the OpenFOAM solver interFoam (color-coded). The physical domain size is 1 m×4 m.

3.1.2 Rising bubble dynamics

We consider a gas bubble rising in a basaltic melt. The bubble, initially at rest, rises due to buoyancy assuming a variety of shapes depending on the system parameters (e.g., liquid viscosity, surface tension, density contrast). Samkhaniani et al. (2012) demonstrated the ability of interFoam to reproduce the different bubble shapes reported in the Grace diagram (Grace1973). Our contribution here focuses on the validation of the solver for bubble stability in magmas through comparison with Suckale et al. (2010a) that used a different numerical method (Fig. 7). In this set of 2D simulations, the main goal is to investigate the ability of the solver to reproduce the breakage of a bubble in relation to the shape that it may assume during its rise. Breakage may occur because of the small, random perturbations that form at the melt–bubble interface. No-slip boundary conditions are used for top and bottom boundaries and periodic conditions for the sides. In the volcanic context, two parameters change significantly with respect to water experiments (Samkhaniani et al.2012): the density ρ and the viscosity μ of the liquid, while for water ρ=103kgm-3 and μ=10-3Pas, even a low-viscosity silicate melt (e.g., basalt) has viscosity values of the order of 10–100 Pa s and the density is above 2500 kg m−3. Surface tension is σ=0.073 N m−1 for water, while a reasonable value for magmas is σ=0.15 N m−1 (Colucci et al.2016). In our tests we set σ=0.3 N m−1 and ρ=3500 kg m−3 to be consistent with Suckale et al. (2010a). With a being the bubble radius and v0 the rise velocity, the relevant nondimensional numbers derived directly from the governing equations for an incompressible melt are as follows (Suckale et al.2010a): Reynolds number Re=ρv0aμ (inertial to viscous forces), Froude number Fr=v0ga (inertia to buoyancy forces), Weber number We=ρv02aσ (inertia to surface tension), and gas to liquid viscosity ratio Π=μg/μ. We can also define the Eötvös number (Eo), which is a combination of Fr and We (Eo=FrWe). Let us note that to be consistent with Suckale et al. (2010a) all nondimensional numbers here are based on the bubble radius instead of the bubble diameter, which is also commonly used in the literature (e.g., Roghair et al.2011). Considering a constant Π=10-6, bubble regimes can be classified using only two adimensional numbers, Re and Eo. The Reynolds and Eötvös numbers control bubble stability, deformation, and breakup. Indicatively, for Eo<1 and Re<1 the bubble is stable and preserves its initial spherical shape even in the presence of small perturbations of its interface. For Eo>1 and Re<1 the bubble deforms and may break up if random perturbations significantly affect its surface, while for Eo>1 and Re>1 breakup occurs invariably.

Figure 7Simulations of bubble rise in a basaltic melt using interFoam are compared with the results of Suckale et al. (2010a) (black lines) for three different regimes: (a) no breakup (Re≈5, Fr2≈0.4, We≈90, and Π=10-6), (b) gradual breakup (Re≈25, Fr2≈0.3, We≈800, and Π=10-6), and (c) catastrophic breakup (Re≈250, Fr2≈0.16, We≈1350, and Π=10-6). For each regime, snapshots at different nondimensional times are shown. To be consistent with Suckale et al. (2010a) here we use the square of the Froude number (Fr2=v02ga).

Overall, breakup mechanisms are well reproduced in our simulations and bubble shapes at given nondimensional times are consistent with those reported by Suckale et al. (2010a) for similar values of the nondimensional numbers (Fig. 7) and similar space resolution (20 cells per bubble radius). For the no-breakup regime (Fig. 7a), the bubble shape in our simulation displays two additional thin wings. In the gradual breakup regime (Fig. 7b) small droplets are formed at the rear of the bubble. The results are reported here with higher resolution (40 cells per bubble radius), since with lower resolution the bubble presents a slightly different shape with a flatter head. In the catastrophic breakup regime (Fig. 7c), the bubble immediately collapses, forming a large number of small- to medium-sized bubbles.

3.2 Diffusive bubble growth

Here we demonstrate the ability of the ODE solver multiComponentODERPShellDStatic to simulate bubble growth in a rhyolitic melt by expansion and mass diffusion. Our solution has been benchmarked by comparison with Lyakhovsky et al. (1996) for the diffusive growth of water gas bubbles under instantaneous decompression of a hydrated rhyolitic melt. To reproduce the reference solution we assumed a quasi-static diffusion in the liquid shell around the bubble. The quasi-static approximation holds when diffusion is fast relative to decompression rate (Lensky et al.2004). The reference solutions for three different values of the diffusion coefficient are well reproduced by our model (Fig. 8). We repeat the same simulations with the addition of 1 % of CO2 (red lines in Fig. 8). The multicomponent saturation surface is calculated using SOLWCAD (Papale et al.2006). The bubble radius increases by ≈50 % and the gas volume fraction triples (see Eq. B6).

Figure 8Temporal evolution of bubble radius for an instantaneous decompression from p0=150 MPa to pL=120 MPa. In blue is the comparison between the MagmaFOAM model multiComponentODERPShellDStatic (solid lines) and numerical solutions from Lyakhovsky et al. (1996) (dashed lines) that practically coincide for three different values of the diffusion coefficient of H2O. The red lines represent the same simulations with 1 wt % of CO2 added in the melt. The diffusion coefficient of CO2 is 1 order of magnitude smaller than H2O. Initial conditions and parameters (see Appendix B) in all simulations are ρL=2300 kg m−3, μ=5×104 Pa s, σ=0.06 N m−1, T=1123 K, pG(t=0)=p0+2σ/R(t=0), R(t=0)=10-7 m, S0=2×10-4 m, and CH2O0=0.053wt%. Saturation concentration is computed using SOLWCAD (Papale et al.2006).

3.3 Eulerian multi-fluid modeling

In this section we test the ability of the OpenFOAM two-fluid Eulerian solver reactingTwoPhaseEulerFoam to deal with flow problems with a large number of small (unresolved) gas bubbles dispersed into a liquid phase. The reactingTwoPhaseEulerFoam solver, coupled with the MagmaFOAM libraries, is tested in problems involving multiphase shock tubes as well as by simulating a multiphase–multicomponent reacting box.

3.3.1 Shock tube simulations

Decompression of a pressurized bubbly magma is a common trigger of explosive volcanic eruptions (e.g., Gonnermann and Manga2007). When a high-pressure magma reservoir is decompressed, a shock wave and a contact wave propagate into the low-pressure region, typically the atmosphere, and a rarefaction wave propagates into the bubbly magma (Koyaguchi and Mitani2005; La Spina et al.2017), akin to shock tube devices. The latter have been extensively used to study wave propagation phenomena in compressible fluids. Usually high- and low-pressure regions are separated by a diaphragm, the instantaneous removal of which initiates highly transient dynamics (Stadtke2006). Assuming strictly one-dimensional flow conditions (i.e., ignoring the effects of shear viscosity), the shock tube mathematically represents a Riemann problem in which the initial velocities on both sides have been set to zero. For the specific case of ideal equation of state, an analytical solution can be derived for a pure single phase (Stadtke2006). Multiphase flow processes are generally governed by deviations from mechanical and thermal equilibrium between the phases. Nevertheless, the assumptions of homogeneous flow (equal phase velocities) and thermal equilibrium between the phases allow us to define a special limiting case for which a semi-analytical solution can be derived (Stadtke2006). We test the reactingTwoPhaseEulerFoam solver in inviscid one-dimensional and viscid axisymmetric simulations of single-phase and two-phase shock tubes. Axisymmetric simulations allow us to investigate the effect of melt viscosity on the radial velocity profile through the Giordano et al. (2008) model. The Lange and Carmichael (1987) equation of state is tested here for the propagation of rarefaction and shock waves.

Single phase. We perform a standard Sod shock tube benchmark for pure air gas using a perfect gas equation of state. Nearly perfect agreement between the simulation and the analytical solution has been found by discretizing the one-dimensional computational domain with 5000 cells. Then, we test the solver by simulating a shock tube with pure liquid water using the SPWAT equation of state (Stadtke2006) implemented in MagmaFOAM. Discretizing the computational domain with the same number of cells, the contact and shock wave discontinuities are well resolved and do not display any instabilities. Finally, we perform a shock tube simulation with pure liquid basalt (Table D1) using the equation of state for silicate melts proposed by Lange and Carmichael (1987) and implemented in MagmaFOAM. We use the same computational domain and the same numerical schemes used in the liquid water test. Across the shock discontinuity the solution is more diffusive compared to the test for liquid water, while the contact discontinuity is still well resolved. The figures for the single-phase shock tubes described above are reported in Appendix C.

Two phases. We perform two-phase shock tube simulations for gas air–liquid water and gas water–liquid basalt (Table D1) shock tubes (Figs. 9, 10, 11). The equations of state for each phase are the same as for the single-phase cases. In all the simulations, the size of the dispersed phase (i.e., gas bubbles or liquid droplets), instead of being determined by a proper model (i.e., bubble growth model), is kept constant and used as a tuning parameter. This unphysical assumption allows us to control the mechanical and thermal disequilibrium between the gas and liquid phases in order to compare the simulation with the limiting analytical solution for a homogeneous flow (Stadtke2006). It is worth noting that, even if the size of the dispersed phase is kept constant, its volume fraction can change according to the mass conservation equations.

First, we benchmark the solver with a gas air–liquid water shock tube for which a limiting analytical solution is provided (Stadtke2006) (Fig. 9). Initial gas volume fraction is 0.1 in the high-pressure region (to the left of the interface) and 0.05 in the low-pressure region (to the right of the interface). Overall, we find remarkably good agreement with the exact solution. The contact and shock wave discontinuities are well resolved and do not display any instabilities. The numerical solution is only slightly diffusive at the onset of the rarefaction wave. The overshoot visible in the velocity is produced by the mechanical decoupling between the liquid and dispersed gas phase (Stadtke2006) and disappears, reducing the bubble size, as will be discussed in the next subsection.

The same simulation setup is used to simulate a basalt (liquid)–water (gas) shock tube (Fig. 10). In this case the simulation is axisymmetric in order to understand the role of melt viscosity. The shape of the axial profiles of pressure, velocity, gas volume fraction, and mixture density are similar to a 1D shock tube (Fig. 9) for the air–water system. Velocity profiles along the radial coordinate are flat with a narrow boundary layer near the walls. In this case the higher viscosity (≈10 Pa s) drastically reduces the mechanical phase decoupling and the phase velocities are superimposed.

Finally, we use the same simulation setup in Figs. 9 and 10, but with an initial gas volume fraction in the low-pressure region (to the right of the interface) equal to 1 (Figs. 11 and C4 in Appendix C). This configuration is more appropriate for a volcanic application in which the shock wave travels in the atmosphere. In this case the continuous and dispersed phases can invert, and thus the bubbly flow, in which bubbles are dispersed in the continuous liquid phase, becomes a particle flow, in which the liquid droplets are dispersed in the gas. This process, usually called fragmentation in volcanological literature, can be modeled as a first approximation using a critical volume fraction criterion (0.5<α<0.7; e.g., Sparks1978, or La Spina et al.2017). When the rarefaction wave propagates into the high-pressure region (i.e., left side), the bubbly magma expands, accelerates, and fragments. Due to the higher compressibility of the gas phase compared to the liquid melt, the temperature subplot shows phase decoupling during expansion, the amount of which depends on the adopted heat transfer model.

Figure 9Results at time t=0.015 s for the air–water shock tube using the SPWAT EOS. Solid lines: OpenFOAM; dashed lines: analytical solution (Stadtke2006); black lines: mixture; blue lines: liquid (water); red lines: gas (air). Mixture density is calculated a posteriori as ρmix=αgρg+(1-αg)ρl, where l is liquid and g is gas. At time 0, the interface dividing the high-pressure (left, l) from low-pressure (right, r) zone is placed at 2.5 m. Initial conditions: Pl=0.5 MPa, Pr=0.1 MPa, and Tl=Tr=300 K for both phases; gas volume fraction: αl=0.05, αr=0.1, and Ul=Ur=0 for both phases. Isobaric heat capacities of gas air and liquid water are Cpg=1004.5 J kg−1 K−1 and Cpl=4195 J kg−1 K−1, respectively (, last access: 27 April 2022). Prandtl numbers of air and water are 0.7 and 2.289, respectively, corresponding to thermal conductivities of 0.02 and 0.67 W K−1 m−1 (, last access: 27 April 2022).


Figure 10Results at time t=0.0065 s for the axisymmetric gas water–basalt shock tube using the Lange–Carmichael EOS (Lange and Carmichael1987) and the viscosity model of Giordano et al. (2008). Blue lines: liquid (basalt); red lines: gas (water). At time 0, the interface dividing the high-pressure (left, l) from low-pressure (right, r) zone is placed at 2.5 m. Initial conditions: Pl=10 MPa, Pr=0.1 MPa, and Tl=Tr=1273 K for both phases; gas volume fraction: αl=0.05, αr=0.1, and Ul=Ur=0 for both phases. Isobaric specific heat capacities in the gas and liquid phase are CPg=2510 J kg−1 K−1 (, last access: 27 April 2022) as well as CPl,H2O=2278Jkg-1K-1 and CPl,basalt=1600Jkg-1K-1, respectively (Lesher and Spera2015). Thermal conductivity of the liquid is 1.5 W K−1 m−1 (Lesher and Spera2015); for the water gas phase a Prandtl number of 0.9 is used, corresponding to a thermal conductivity of 0.14 W K−1 m−1 (, last access: 27 April 2022).

Figure 11Results at time t=0.001 s for the axisymmetric gas water–basalt shock tube using the Lange–Carmichael EOS (Lange and Carmichael1987) and the viscosity model of Giordano et al. (2008), with liquid phase switching from continuous to dispersed. Blue lines: liquid (basalt); red lines: gas (water). At time 0, the interface dividing the high-pressure (left, l) from low-pressure (right, r) zone is placed at 2.5 m. Initial conditions: Pl=10 MPa, Pr=0.1 MPa, and Tl=Tr=1273 K for both phases; gas volume fraction: αl=0.05, αr=1, and Ul=Ur=0 for both phases. Isobaric specific heat capacities in the gas and liquid phase are CPg=2510Jkg-1K-1 (, last access: 27 April 2022) as well as CPl,H2O=2278Jkg-1K-1 and CPl,basalt=1600Jkg-1K-1, respectively (Lesher and Spera2015). Thermal conductivity of the liquid is 1.5WK-1m-1 (Lesher and Spera2015); for the water gas phase a Prandtl number of 0.9 is used, corresponding to thermal conductivity of 0.14WK-1m-1 (, last access: 27 April 2022).

The phase coupling problem. Even if we limit to bubbly flow regimes, magmatic systems are characterized by a wide range of viscosities (from 0.1 to 109Pa s) and bubble sizes (from a few microns to decimeters). The bubble relaxation time (τb) is directly proportional to the square of the bubble diameter and inversely proportional to the kinematic viscosity of the continuous liquid phase (τbDb2/νl). In magmatic phenomena, when considering small bubbles (e.g., 100 µm) and even relatively low viscosities (e.g., 10 Pa s), τb can reach very small values (τb10-6s), resulting in very strong mechanical phase coupling. Numerical algorithms like the one implemented in OpenFOAM, based on segregated solvers, require special care in order to ensure convergence of the solution when phase coupling is tight (Karema and Lo1999). The Partial Elimination Algorithm (PEA), implemented in OpenFOAM, shows the best convergence performance compared to other algorithms (Karema and Lo1999; Venier et al.2016). Here we test the PEA method for shock tube simulation conditions within the range of interest for magmatic applications. In Fig. 12 we compare the analytical solution to the simulation results for different values of τb obtained by changing the bubble diameter. Decreasing τb from 10−1 to 10−3 s, the velocities of the two phases tend to overlap as expected and agree with the homogeneous analytical solution. However, by further decreasing τb the solution diverges even when increasing the number of corrector cycles 40 times. This is an important limitation in the use of the multi-fluid solver. A possible workaround, currently under investigation, is to implement a limiter for the relaxation time.

Figure 12Air–water shock tube simulations at different relaxation times τ and number of correctors. Dashed lines: analytical solution; solid lines: simulation. Blue lines: liquid (water); red lines: gas (air). At time 0, the interface dividing the high-pressure (left, l) from low-pressure (right, r) zone is placed at 2.5 m. Initial conditions: Pl=0.5 MPa, Pr=0.1 MPa, and Tl=Tr=300K for both phases; gas volume fraction: αl=0.05, αr=0.1, and Ul=Ur=0 for both phases. Isobaric heat capacities of gas air and liquid water are CPg=1004.5Jkg-1K-1 and CPl=4195Jkg-1K-1, respectively (, last access: 27 April 2022). Prandtl numbers of air and water are 0.7 and 2.289, respectively, corresponding to thermal conductivities of 0.02 and 0.67 W K−1 m−1, respectively (, last access: 27 April 2022).

3.3.2 Reacting box

A many-bubble system at zero gravity, in which bubbles grow by mass diffusion, is analyzed here. The liquid phase is a basaltic melt (Table D1) with dissolved water and carbon dioxide whose properties are modeled by the Lange and Carmichael (1987) equation of state and the rheological equation of Giordano et al. (2008). The ideal gas equation of state has been used for the gas phase. As this is a many-bubble system, bubble growth is approximated through a subgrid model (see Sect. 2.2). The mass transfer coefficient (i.e., ki in Eq. 4) is calculated according to a spherical model and the heat transfer coefficient according to spherical and Ranz–Marshall models, which are both already implemented in OpenFOAM. In addition, the one-group interfacial area transport equation (IATE) model (Ishi and Hibiki2006) is used to compute the interfacial area required by the mass and heat transfer rate. The IATE is a fundamental equation, formulated from the Boltzmann transport equation, describing the change in surface area between phases, assuming a spherical shape of the dispersed phase (Ishi and Hibiki2006).

Figure 13Reacting box simulation. At time 0 a small amount of gas (volume fraction αg=10-4) is uniformly distributed in the box, and the basaltic melt is oversaturated in H2O (5 wt %) and CO2 (5000 ppm) at 80 MPa and 1373 K. The diffusion coefficients for H2O and CO2 in the basalt are 10−9 and 10-10m2s-1, respectively (Baker et al.2005). Isobaric specific heat capacities in the gas phase are CPg,H2O=2900Jkg-1K-1 and CPg,CO2=1390Jkg-1K-1 (, last access: 27 April 2022, Beaton et al.1987) and in the liquid phase CPl,H2O=2278Jkg-1K-1, CPl,CO2=1600Jkg-1K-1, and CPbasalt=1600Jkg-1K-1 (Lesher and Spera2015). Thermal conductivity of the liquid is 1.5WK-1m-1 (Lesher and Spera2015); for the gas phase Prandtl numbers of 0.9 for H2O and 0.7 for CO2 are used, corresponding to thermal conductivities of 0.16 and 0.09 WK-1m-1 (, last access: 27 April 2022, Beaton et al.1987).

At time zero, a small amount of gas is uniformly distributed in the box and the liquid–gas system is out of thermodynamic equilibrium because the liquid is oversaturated in both H2O and CO2. After 4×105s, H2O has reached thermodynamic equilibrium, increasing the gas volume fraction to ≈33 % (Fig. 13) and the bubble size from 1 to about 4 cm (not shown in the figure). This time is consistent with the timescale that characterizes diffusive mass transfer of H2O (diffusion coefficient D=10-9m2s-1, Baker et al.2005) around a bubble with radius R≈2 cm (τd=R2/D; Lensky et al.2004). The dissolved CO2 is still out of thermodynamic equilibrium, as expected, because the diffusion coefficient is lower, being set to D=10-10m2s-1 (Baker et al.2005). The density and viscosity of the liquid vary with the decreasing dissolved H2O. The gas density decreases because of increasing H2O and decreasing CO2 that produce a decrease in the molar mass of the gas mixture.

4 Conclusions

In this work we have presented MagmaFOAM, a library based on OpenFOAM that contains dedicated tools for the simulation of multiphase flows in magmatic systems. The MagmaFOAM implementation results in a flexible framework which is ideal for development, testing, coupling, and application of the great collection of existing and future modeling strategies needed to tackle the variety of nonlinear multiscale processes characterizing magma flows. MagmaFOAM includes dedicated multicomponent constitutive models for dealing with realistic properties for silicate melt–gas systems as well as different utilities that automatize pre- and post-processing operations. We have analyzed a number of test cases that illustrate the current capabilities and limitations of different modeling approaches in solving well-defined and reproducible flow problems, establishing solid ground for future model selection, improvement, and intercomparison studies. We have shown some of the ingredients that can be used for simulating the interaction among different silicate melts, as well as between melts and fluid phases, under different assumptions and aimed at different portions of the magmatic system (deep reservoirs vs. shallow conduits). Applications of MagmaFOAM can thus include the study of magma mingling and mixing, as well as slug rising dynamics or volatile flushing. Nevertheless, important limitations remain, most notably the development of a magma-specific mixture approach or the intrinsic complications in modeling the transition from tight to loose phase coupling (Sect. 3.3.1).

The framework described in this work allows for maximum flexibility and adaptability, giving the possibility to explore magmatic systems with different approaches given the specific conditions aimed at. As an example, the MagmaFOAM modular approach allows the coupling of its bubble growth models with both single- and multi-fluid solvers, Lagrangian particle tracking, or with more complex constitutive equations. Indeed, at different stages within the evolution of magmatic plumbing systems, different modeling approaches can be more appropriate to capture the fundamental physics governing the dynamics: while low-gas-fraction, deep reservoirs may well be approximated by mixture theory, at shallower levels phase decoupling becomes important and multi-fluid descriptions are more appropriate.

The tool is meant to be under continuous development and is already underway. The addition of population balance equations to single- and multi-fluid models to statistically describe the dispersed phases (bubbles and crystals, Marchisio and Fox2013) will improve our understanding of how polydispersity can impact magmatic system evolution (Colucci et al.2017a; de' Michieli Vitturi and Pardini2021). In large-scale multi-fluid simulations, the exchanges of mass, momentum, and energy through the interface between phases need to be modeled accurately to determine the rate of phase change and the degree of mechanical and thermal disequilibrium between phases. The population balance is a statistical approach for modeling the mesoscale dynamics widely used in chemical engineering, which describe the temporal and space evolution of a large number of particles through a number distribution function (Yeoh and Tu2019a). In this way microscopic processes involving bubble dynamics and interactions between bubbles can be included in large-scale multi-fluid simulations. In fact, DNS allows modeling particle–particle interactions and capturing emerging behaviors in complex systems; however, the large quantity of microphysics taken into account in DNS has to be filtered and condensed in a sub-model to be used in large-scale simulations. Mesoscopic models represent intermediate models that describe, through a set of mesoscale variables, the microphysics of the system. The formulation of population balance requires adequate closure models for the microphysics that can be developed with the aid of experimental (Mancini et al.2016) and DNS investigations (Marchisio and Fox2013). The inclusion of Lagrangian tracers will result in a more detailed description, with respect to multi-fluid models, of the microphysics that determine the macroscopic properties driving the dynamics. In the Eulerian–Lagrangian approach, bubbles are treated as discrete Lagrangian particles in an ambient Eulerian continuous flow (e.g., Ghahramani et al.2019). This approach in fact is more appropriate than multi-fluid models when the number of particles is too small to be treated as a continuum or when the behavior of single particles (e.g., rapidly expanding or contracting bubbles) is so specific that they are not well represented by unique averaged field density, velocity, or temperature (e.g., Ghahramani et al.2019). With respect to the DNS approach, whereby bubble–bubble and bubble–melt interactions emerge self-consistently, in Eulerian–Lagrangian models phase interactions are defined by constitutive models. However, the Eulerian–Lagrangian approach, compared to the DNS, allows simulating larger populations of particles at a much lower computational cost. The study of complex mixing behavior in magma mushes is only an example of possible applications (Bergantz et al.2015). Finally, engineering applications have benefited from models that combine different approaches, e.g., interface-resolving and subgrid dispersed phase modeling with single- or multi-fluid frameworks. These hybrid models, although not fully mature yet, in principle allow modeling the broad range of interface scales that typically characterize gas–liquid flows including regime transitions at the same time (e.g., Wardle and Weller2013). From a volcanological perspective, predicting flow regime changes is of crucial importance to understand effusive–explosive transitions in eruptive activity (Gonnermann and Manga2007).

Appendix A: Linear Rayleigh–Taylor instability

Figure A1b shows how a small wave number perturbation (k=0.15) initially grows with a slower nonconstant growth rate. Overall this effect makes the extrapolated growth rate smaller than expected. However, after a relatively small time interval, the growth rate becomes constant with a value that results to be in good agreement with the theoretical one (Fig. A2). This spurious effect gradually decreases until it disappears as the wave number of the perturbation increases (Fig. A1a). The simulations are done using the solver interFoam with adaptable time step (Comax=0.01).

Figure A1Time evolution of the amplitude of two single-mode perturbations (k=0.5a; k=0.15b) for the linear Rayleigh–Taylor instability benchmark. The growth rate of the perturbation is extrapolated with a linear regression excluding data in late (physical) and eventually early (spurious) phases characterized by nonlinear effects (data not marked with an asterisk).


Figure A2Extrapolated growth rate for two perturbations with linear regression excluding (blue) or not excluding (red) data in the initial phase characterized by nonlinear spurious effects.


Appendix B: multiComponentODERPShellDStatic: model equations

A modified form of the Rayleigh–Plesset equation describes the hydrodynamics of the growth of a multicomponent spherical bubble in a finite incompressible shell of liquid.


In the above, ρL is the liquid density, R is the bubble radius, S is the radius of the shell (S3=S03+R3; S03=S3(t=0)-R3(t=0)), pG is the gas pressure inside the bubble, pL is the pressure acting on the outside of the liquid shell, σ is the surface tension, t is the time, and μ is the liquid dynamic viscosity that depends on the concentration of dissolved volatiles in the shell. Given pL(t) this represents an equation that can be solved to find R(t) provided pG(t). pG is given by combining the mass conservation of the gas phase with an equation of state for a perfect gas. Mass conservation of the gas phase is given by

(B2) d d t R 3 ρ G = 3 R 2 ρ L i = 1 N D i C i r r = R ,

where Di is the mass diffusivity and Ci the concentration of the ith species dissolved in the melt. Mass conservation of the ith dissolved species is given by

(B3) d d t R 3 ρ G Y i = 3 R 2 ρ L D i C i r r = R ,

where Yi is the concentration in the gas phase. Assuming local thermodynamic equilibrium at the bubble–melt interface (i.e., Ci(R)=Cisat(pG), where Cisat is the saturation concentration), zero gradient boundary conditions at the shell boundary, and a quasi-static diffusion in the shell (Lyakhovsky et al.1996), the term in square brackets in Eqs. (B2) and (B3) is given by

(B4) d C i d r r = R = S 0 3 C i 0 - C i sat + 1 ρ L R 3 ρ G Y i t = 0 - R 3 ρ G Y i S 0 3 R - 3 2 S 2 - R 2 R 2 ,

where Ci0 is the concentration in the melt at time 0. Assuming constant viscosity, the term in Eq. (B1) is analytically integrated to obtain

(B5) - 3 R S μ ( C i ( r ) ) r 4 d r = μ 1 S 3 - 1 R 3 .

For a monodispersed distribution, the gas volume fraction is given by

(B6) α = R 3 S 3 .
Appendix C: Shock Tube

Figures C1, C2, and C3 show results from the single-phase shock tube simulations discussed in Sect. 3.3.1. Figure C4 shows results from the air–water shock tube with liquid phase switching from continuous to dispersed.

Figure C1Results at time t=0.006 s for the air Sod shock tube. Dashed lines: analytical solution; solid lines: simulation. At time 0, the interface dividing the high-pressure (left, l) from low-pressure (right, r) zone is placed at 0 m. Initial conditions: Pl=0.1 MPa, Pr=0.01 MPa; Tl=348.432 K, Tr=278.746 K; gas volume fraction αl=1, αr=0; Ul=Ur=0. Isobaric heat capacity is CP=1004.5Jkg-1K-1, corresponding to heat capacity ratio γ=1.4.


Figure C2Results at time t=0.00164 s for single-phase shock tube with liquid water using SPWAT EOS. At time 0, the interface dividing the high-pressure (left, l) from low-pressure (right, r) zone is placed at 0. Initial conditions: Pl=10 MPa, Pr=0.1 MPa; Tl=Tr=300K; gas volume fraction αl=αr=0; Ul=Ur=0. Isobaric heat capacity is CP=4195Jkg-1K-1 (, last access: 27 April 2022).

Figure C3Results at time t=0.002 s for single-phase shock tube with basaltic melt using Lange–Carmichael EOS. At time 0, the interface dividing the high-pressure (left, l) from low-pressure (right, r) zone is placed at 0. Initial conditions: Pl=10 MPa, Pr=0.1 MPa; Tl=Tr=1373K; gas volume fraction αl=αr=0; Ul=Ur=0. Isobaric heat capacity is CP=1600Jkg-1K-1 (Lesher and Spera2015). Thermal conductivity is 1.5WK-1m-1 (Lesher and Spera2015).

Figure C4Results at time t=0.0015 s for the air–water shock tube with liquid phase switching from continuous to dispersed. Blue lines: liquid (water); red lines: gas (air). At time 0, the interface dividing the high-pressure (left, l) from low-pressure (right, r) zone is placed at 2.5 m. Initial conditions: Pl=0.5 MPa, Pr=0.1 MPa; Tl=Tr=300K for both phases; gas volume fraction αl=0.05, αr=1; Ul=Ur=0 for both phases. Isobaric heat capacities of gas air and liquid water are CPg=1004.5Jkg-1K-1 and CPl=4195Jkg-1K-1, respectively (, last access: 27 April 2022). Prandtl numbers of air and water are 0.7 and 2.289, respectively, corresponding to thermal conductivities of 0.02 and 0.67 WK-1m-1 (, last access: 27 April 2022).

Appendix D: Magmatic compositions

Table D1 reports the compositions in terms of major oxides of the magmas used in the simulations shown in the paper.

Table D1Oxide compositions for the magmas used in benchmarking simulations. Amounts are relative.

Download Print Version | Download XLSX

Code and data availability

The version of the model used to produce the results shown in this paper, as well as input data and scripts to replicate all the simulations presented in this paper, is archived on Zenodo ( et al.2021).

Author contributions

FB and SC developed and tested the software including pre- and post-processing, performed the simulations, and wrote the first draft of the paper. CPM contributed to paper writing and supervised the work. JM worked on bubble dynamics and their analysis. MdMV provided guidance and help on the multi-fluid solver. PP provided supervision and reviewed the original draft.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors thank Jenny Suckale and an anonymous reviewer for their constructive comments that considerably improved the paper.

Financial support

This research has been supported by the Istituto Nazionale di Geofisica e Vulcanologia (INGV) and Istituto Nazionale di Oceanografia e Geofisica Sperimentale (OGS) under the HPC-TRES program award (grant no. 2016-05) to Federico Brogi. This research has received funding from European Union's Horizon 2020 research and innovation program under the EUROVOLC project (grant no. 731070) and under the ChEESE project (grant no. 823844) as well as Italy's MIUR PRIN (grant no. 2015L33WAK) and from INGV Pianeta Dinamico (CHOPIN).

Review statement

This paper was edited by Lutz Gross and reviewed by Jenny Suckale and one anonymous referee.


Annen, C., Blundy, J. D., and Sparks, R. S. J.: The Genesis of Intermediate and Silicic Magmas in Deep Crustal Hot Zones, J. Petrol., 47, 505–539,, 2006. a

Aulisa, E., Manservisi, S., and Scardovelli, R.: A mixed markers and volume-of-fluid method for the reconstruction and advection of interfaces in two-phase and free-boundary flows, J. Computat. Phys., 188, 611–639, 2003. a

Bachmann, O. and Bergantz, G. W.: Rejuvenation of the Fish Canyon magma body: A window into the evolution of large-volume silicic magma systems, Geology, 31, 789–792, 2003. a

Bagagli, M., Montagna, C. P., Papale, P., and Longo, A.: Signature of magmatic processes in strainmeter records at Campi Flegrei (Italy), Geophys. Res. Lett., 44, 718–725,, 2017. a, b

Baker, D. R., Freda, C., Brooker, R. A., and Scarlato, P.: Volatile diffusion in silicate melts and its effects on melt inclusions, Ann. Geophys., 48, 699–717, 2005. a, b, c

Beaton, C., Edwards, D., and Schlünder, E.: Heat Exchanger Design Handbook: HEDH, Physical properties, edited by: Faulkner, L. L., v. 5, VDI-Verlag, ISBN 0-8247-9787-6, 1987. a, b

Bergantz, G., Schleicher, J., and Burgisser, A.: Open-system dynamics and mixing in magma mushes, Nat. Geosci., 8, 793,, 2015. a, b

Brennen, C. E. and Brennen, C. E.: Fundamentals of multiphase flow, Cambridge University Press, ISBN 9780521848046, 2005. a

Brogi, F., Colucci, S., and Montagna, C. P.: MagmaFOAM-1.0: a modular framework for the simulation of magmatic systems, Zenodo [code],, 2021. a

Burgisser, A., Alletti, M., and Scaillet, B.: Simulating the behavior of volatiles belonging to the C–O–H–S system in silicate melts under magmatic conditions with the software D-Compress, Comput. Geosci., 79, 1–14, 2015. a

Burnham, C. W. and Davis, N.: The role of H2O in silicate melts; II, Thermodynamic and phase relations in the system NaAlSi3O8-H2O to 10 kilobars, 700 degrees to 1100 degrees C, Am. J. Sci., 274, 902–940, 1974. a

Cerminara, M., Esposti Ongaro, T., and Berselli, L. C.: ASHEE-1.0: a compressible, equilibrium–Eulerian model for volcanic ash plumes, Geosci. Model Dev., 9, 697–730,, 2016. a

Chandrasekhar, S.: The character of the equilibrium of an incompressible heavy viscous fluid of variable density, in: Mathematical Proceedings of the Cambridge Philosophical Society, Cambridge University Press, 51, 162–178,, 1955. a

Chandrasekhar, S.: Hydrodynamic and hydromagnetic stability, Courier Corporation, in: Dover Books on Physics Series, Dover Publications, ISBN 9780486640716, 2013. a

Cheng, L., Costa, F., and Bergantz, G.: Linking fluid dynamics and olivine crystal scale zoning during simulated magma intrusion, Contrib. Mineral. Petr., 175, 1–14, 2020. a

Chouet, B., Dawson, P., and Nakano, M.: Dynamics of diffusive bubble growth and pressure recovery in a bubbly rhyolitic melt embedded in an elastic solid, J. Geophys. Res.-Sol. Ea., 111,, 2006. a

Colucci, S., Battaglia, M., and Trigila, R.: A thermodynamical model for the surface tension of silicate melts in contact with H2O gas, Geochim. Cosmochim. Ac., 175, 127,, 2016. a

Colucci, S., de' Michieli Vitturi, M., and Landi, P.: CrystalMom: a new model for the evolution of crystal size distributions in magmas with the quadrature-based method of moments, Contrib. Mineral. Petrol., 172, 100,, 2017a. a

Colucci, S., Papale, P., and Montagna, C. P.: Non-Newtonian flow of bubbly magma in volcanic conduits, J. Geophys. Res.-Sol. Ea., 122, 1789–1804,, 2017b. a

Costa, A., Suzuki, Y. J., Cerminara, M., Devenish, B. J., Esposti Ongaro, T., Herzog, M., Van Eaton, A. R., Denby, L. C., Bursik, M., de' Michieli Vitturi, M., Engwell, S., Neri, A., Barsotti, S., Folch, A., Macedonio, G., Girault, F., Carazzo, G., Tait, S., Kaminski, E., Mastin, L. G., Woodhouse, M. J., Phillips, J. C., Hogg, A. J., Degruyter, W., and Bonadonna, C.: Results of the eruptive column model inter-comparison study, J. Volcanol. Geoth. Res., 326, 2–25, 2016. a

Coumans, J., Llewellin, E., Wadsworth, F., Humphreys, M., Mathias, S., Yelverton, B., and Gardner, J.: An experimentally validated numerical model for bubble growth in magma, J. Volcanol. Geoth. Res., 402, 107002,, 2020. a, b

Cussler, E. L. (Ed.): Diffusion mass transfer in fluid systems, in: Cambridge Series in Chemical Engineering, 3rd edn., Cambridge University Press,, 2009. a

de' Michieli Vitturi, M. and Pardini, F.: PLUME-MoM-TSM 1.0.0: a volcanic column and umbrella cloud spreading model, Geosci. Model Dev., 14, 1345–1377,, 2021. a

de' Michieli Vitturi, M., Clarke, A., Neri, A., and Voight, B.: Effects of conduit geometry on magma ascent dynamics in dome-forming eruptions, Earth Planet. Sc. Lett., 272, 567–578, 2008a. a

de' Michieli Vitturi, M., Clarke, A., Neri, A., and Voight, B.: Effects of conduit geometry on magma ascent dynamics in dome-forming eruptions, Earth Planet. Sc. Lett., 272, 567–578,, 2008b. a

de' Michieli Vitturi, M., Neri, A., and Barsotti, S.: PLUME-MoM 1.0: A new integral model of volcanic plumes based on the method of moments, Geosci. Model Dev., 8, 2447–2463,, 2015. a

Deshpande, S. S., Anumolu, L., and Trujillo, M. F.: Evaluating the performance of the two-phase flow solver interFoam, Computat. Sci. Discov., 5, 014016,, 2012. a, b

Druitt, T. H., Costa, F., Deloule, E., Dungan, M., and Scaillet, B.: Decadal to monthly timescales of magma transfer and reservoir growth at a caldera volcano, Nature, 482, 77–80,, 2012. a

Dufek, J.: The Fluid Mechanics of Pyroclastic Density Currents, Annu. Rev. Fluid Mech., 48, 459–485,, 2016. a, b

Dufek, J. and Bergantz, G.: Transient two-dimensional dynamics in the upper conduit of a rhyolitic eruption: A comparison of closure models for the granular stress, J. Volcanol. Geoth. Res., 143, 113–132, 2005. a

Duff, R., Harlow, F., and Hirt, C.: Effects of diffusion on interface instability between gases, Phys. Fluids, 5, 417–425, 1962. a

Esposti Ongaro, T., Cavazzoni, C., Erbacci, G., Neri, A., and Salvetti, M.: A parallel multiphase flow code for the 3D simulation of explosive volcanic eruptions, Parallel Computing, 33, 541–560,, 2007. a

Fang, J., Cambareri, J. J., Rasquin, M., Gouws, A., Balakrishnan, R., Jansen, K. E., and Bolotnov, I. A.: Interface tracking investigation of geometric effects on the bubbly flow in PWR subchannels, Nucl. Sci. Eng., 193, 46–62, 2019. a

Ferry, J. and Balachandar, S.: A fast Eulerian method for disperse two-phase flow, Int. J. Multiphas. Flow, 27, 1199–1226, 2001. a

Garg, D., Papale, P., Colucci, S., and Longo, A.: Long-lived compositional heterogeneities in magma chambers, and implications for volcanic hazard, Sci. Rep.-UK, 9, 3321,, 2019. a, b, c, d, e

Ghahramani, E., Arabnejad, M. H., and Bensow, R. E.: A comparative study between numerical methods in simulation of cavitating bubbles, Int. J. Multiphas. Flow, 111, 339–359,, 2019. a, b

Giordano, D., Russell, J. K., and Dingwell, D. B.: Viscosity of magmatic liquids: a model, Earth Planet. Sc. Lett., 271, 123–134, 2008. a, b, c, d, e, f

Gonnermann, H. M. and Manga, M.: The Fluid Mechanics Inside a Volcano, Annu. Rev. Fluid Mech., 39, 321–356,, 2007. a, b, c

Grace, J.: Shapes and velocities of bubbles rising in infinite liquid, T. I. Chem. Eng.-Lond., 51, 116–120, 1973. a

Griffiths, R. W.: The Dynamics of Lava Flows, Annu. Rev. Fluid Mech., 32, 477–518,, 2000. a, b

He, X., Chen, S., and Zhang, R.: A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh–Taylor instability, J. Computat. Phys., 152, 642–663, 1999. a, b, c, d

Huber, C., Su, Y., Nguyen, C., Parmigiani, A., Gonnermann, H. M., and Dufek, J.: A new bubble dynamics model to study bubble growth, deformation, and coalescence, J. Geophys. Res.-Sol. Ea., 119, 216–239, 2014. a, b

Ishi, M. and Hibiki, T.: Thermo-fluid dynamics of two-phase flow, 2006. a, b, c

Issa, R. I.: Solution of the implicitly discretised fluid flow equations by operator-splitting, J. Computat. Phys., 62, 40–65, 1986. a

Jackson, M. D., Blundy, J. D., and Sparks, R. S. J.: Chemical differentiation, cold storage and remobilization of magma in the Earth's crust, Nature, 564, 405–409,, 2018. a

Jellinek, a. M., Kerr, R. C., and Griffiths, R. W.: Mixing and compositional stratification produced by natural convection: 1. Experiments and their application to Earth's core and mantle, J. Geophys. Res., 104, 7183,, 1999. a

Karema, H. and Lo, S.: Efficiency of interphase coupling algorithms in fluidized bed conditions, Comput. Fluids, 28, 323–360,, 1999. a, b

Keller, T. and Suckale, J.: A continuum model of multi-phase reactive transport in igneous systems, Geophys. J. Int., 219, 185–222,, 2019. a, b, c

Koyaguchi, T. and Mitani, N. K.: A theoretical model for fragmentation of viscous bubbly magmas in shock tubes, J. Geophys. Res.-Sol. Ea., 110, B10202,, 2005. a

La Spina, G., de' Michieli Vitturi, M., and Clarke, A.: Transient numerical model of magma ascent dynamics: application to the explosive eruptions at the Soufrière Hills Volcano, J. Volcanol. Geoth. Res., 336, 118–139, 2017. a, b

Lange, R. A. and Carmichael, I. S.: Densities of Na2O-K2O-CaO-MgO-FeO-Fe2O3-Al2O3-TiO2-SiO2 liquids: new measurements and derived partial molar properties, Geochim. Cosmochim. Ac., 51, 2931–2946, 1987. a, b, c, d, e, f, g, h

Lensky, N. G., Lyakhovsky, V., and Navon, O.: Radial variations of melt viscosity around growing bubbles and gas overpressure in vesiculating magmas, Earth Planet. Sc. Lett., 186, 1–6, 2001. a

Lensky, N. G., , Navon, O., and Lyakhovsky, V.: Bubble growth during decompression of magma: experimental and theoretical investigation, J. Volcanol. Geoth. Res., 129, 7–22, 2004. a, b, c

Lesher, C. E. and Spera, F. J.: Thermodynamic and transport properties of silicate melts and magma, in: The Encyclopedia of Volcanoes (Second Edition), edited by: Sigurdsson, H., Elsevier, 113–141,, 2015. a, b, c, d, e, f, g, h, i, j, k

Longo, A., Papale, P., Vassalli, M., Saccorotti, G., Montagna, C. P., Cassioli, A., Giudice, S., and Boschi, E.: Magma convection and mixing dynamics as a source of ultra-long-period oscillations, Bull. Volcanol., 74, 873–880, 2012. a

Lyakhovsky, V., Hurwitz, S., and Navon, O.: Bubble growth in rhyolitic melts: experimental and numerical investigation, Bull. Volcanol., 58, 19–32,, 1996. a, b, c, d

Macedonio, G., Neri, A., Martí, J., and Folch, A.: Temporal evolution of flow conditions in sustained magmatic explosive eruptions, J. Volcanol. Geoth. Res., 143, 153–172,, 2005. a

Mancini, S., Forestier-Coste, L., Burgisser, A., James, F., and Castro, J.: An expansion–coalescence model to track gas bubble populations in magmas, J. Volcanol. Geoth. Res., 313, 44–58, 2016. a

Manga, M.: Waves of bubbles in basaltic magmas and lavas, J. Geophys. Res.-Sol. Ea., 101, 17457–17465, 1996. a

Marchisio, D. L. and Fox, R. O. (Eds.): Multiphase reacting flows: modelling and simulation, Springer,, 2007. a

Marchisio, D. L. and Fox, R. O.: Computational Models for Polydisperse Particulate and Multiphase Systems, Cambridge Series in Chemical Engineering, Cambridge University Press,, 2013. a, b

Marschall, H.: Towards the numerical simulation of multi-scale two-phase flows, PhD thesis, Technische Universität München,, 2011. a, b

Martí, J., Zafrilla, S., Andújar, J., Jiménez-Mejías, M., Scaillet, B., Pedrazzi, D., Doronzo, D., and Scaillet, S.: Controls of magma chamber zonation on eruption dynamics and deposits stratigraphy: The case of El Palomar fallout succession (Tenerife, Canary Islands), J. Volcanol. Geoth. Res., 399, 106908,, 2020. a

Massol, H. and Koyaguchi, T.: The effect of magma flow on nucleation of gas bubbles in a volcanic conduit, J. Volcanol. Geoth. Res., 143, 69–88,, 2005. a

Mauro, J. C., Yue, Y., Ellison, A. J., Gupta, P. K., and Allan, D. C.: Viscosity of glass-forming liquids, P. Natl. Acad. Sci. USA, 106, 19780–19784,, 2009. a

Melnik, O.: Dynamics of two-phase conduit flow of high-viscosity gas-saturated magma: large variations of sustained explosive eruption intensity, B. Volcanol., 62, 153–170,, 2000. a

Melnik, O. and Sparks, S.: Transient models of conduit flows during volcanic eruptions, Statistics in Volcanology, pp. 201–214, 2006. a

Mitchner, M. and Landshoff, R.: Rayleigh-Taylor Instability for Compressible Fluids, Phys. Fluids, 7, 862–866, 1964. a

Moin, P. and Mahesh, K.: Direct numerical simulation: a tool in turbulence research, Annu. Rev. Fluid Mech., 30, 539–578, 1998. a

Montagna, C., Papale, P., and Longo, A.: Timescales of mingling in shallow magmatic reservoirs, Geological Society, London, Special Publications, 422, SP422–6, 2015. a, b

Moreno Soto, Ã., Enríquez, O. R., Prosperetti, A., Lohse, D., and van der Meer, D.: Transition to convection in single bubble diffusive growth, J. Fluid Mech., 871, 332–349,, 2019. a

Morgavi, D., Arienzo, I., Montagna, C. P., Perugini, D., and Dingwell, D. B.: Magma Mixing: History and Dynamics of an Eruption Trigger, in: Volcanic Unrest, 123–137,, 2017. a

Neri, A., Esposti Ongaro, T., Macedonio, G., and Gidaspow, D.: Multiparticle simulation of collapsing volcanic columns and pyroclastic flow, J. Geophys. Res.-Sol. Ea., 108, B4, 2003. a

Newman, S. and Lowenstern, J. B.: VolatileCalc: a silicate melt–H2O–CO2 solution model written in Visual Basic for excel, Comput. Geosci., 28, 597–604,, 2002. a

Ongaro, T. E., Cavazzoni, C., Erbacci, G., Neri, A., and Salvetti, M.-V.: A parallel multiphase flow code for the 3D simulation of explosive volcanic eruptions, Parallel Computing, 33, 541–560, 2007. a

Papale, P.: Modeling of the solubility of a two-component H2O+ CO2 fluid in silicate liquids, American Mineralogist, 84, 477–492, 1999. a

Papale, P.: Dynamics of magma flow in volcanic conduits with variable fragmentation efficiency and nonequilibrium pumice degassing, J. Geophys. Res., 106, 11043,, 2001. a, b

Papale, P., Moretti, R., and D., B.: The compositional dependence of the saturation surface of H2O+CO2 fluids in silicate melts, Chem. Geol., 229, 78–95, 2006. a, b, c, d, e, f

Parmigiani, A., Huber, C., and Bachmann, O.: Mush microphysics and the reactivation of crystal-rich magma reservoirs, J. Geophys. Res.-Sol. Ea., 119, 6308–6322, 2014. a

Perugini, D. and Poli, G.: The mixing of magmas in plutonic and volcanic environments: Analogies and differences, Lithos, 153, 261–277,, 2012. a

Perugini, D., Poli, G., Petrelli, M., De Campos, C., and Dingwell, D.: Time-scales of recent Phlegrean Fields eruptions inferred from the application of a “diffusive fractionation” model of trace elements, Bul. Volcanol., 72, 431–447,, 2010. a

Proussevitch, A. A. and Sahagian, D. L.: Dynamics and energetics of bubble growth in magmas: Analytical formulation and numerical modeling, J. Geophys. Res., 103, 18223–18251, 1998. a, b

Proussevitch, A. A., Sahagian, D. L., and Anderson, A. T.: Dynamics of Diffusive Bubble Growth in Magmas: Isothermal Case, J. Geophys. Res., 98, 22283–22307, 1993. a

Qin, Z. and Suckale, J.: Flow-to-Sliding Transition in Crystal-Bearing Magma, J. Geophys. Res.-Sol. Ea., 125, e2019JB018549,, 2020. a, b

Rider, W. J. and Kothe, D. B.: Reconstructing volume tracking, J. Computat. Phys., 141, 112–152, 1998. a

Roenby, J., Larsen, B. E., Bredmose, H., and Jasak, H.: A new volume-of-fluid method in openfoam, in: VII International Conference on Computational Methods in Marine Engineering, Nantes: International Center for Numerical Methods in Engineering, Proceedings of the VI International Conference on Computational Methods in Marine Engineering, 15–17 May 2017, Nantes, France, 266–277, ISBN 978-84-943928-6-3, 2017. a

Roghair, I., Lau, Y., Deen, N., Slagter, H., Baltussen, M., Van Sint Annaland, M., and Kuipers, J.: On the drag force of bubbles in bubble swarms at intermediate and high Reynolds numbers, Chemical Engineering Science, 66, 3204–3211,, 2011. a

Ruprecht, P., Bergantz, G. W., and Dufek, J.: Modeling of gas-driven magmatic overturn: Tracking of phenocryst dispersal and gathering during magma mixing, Geochem. Geophys. Geosyst., 9, 7,, 2008. a, b, c, d

Sahagian, D. and Proussevitch, A.: Standardized model runs and sensitivity analysis using the “Bubbledrive-1” volcanic conduit flow model, J. Volcanol. Geoth. Res., 143, 173–185,, 2005. a

Samkhaniani, N., Ajami, A., Kayhani, M. H., and Dari, A. S.: Direct numerical simulation of single bubble rising in viscous stagnant liquid, in: International Conference on Merchanical, Automobile and Robotics Engineering, International Conference on Mechanical, Automobile and Robotics Engineering (ICMAR), January 2012, Penang, Malaysia, (last access: 4 May 2022), 2012. a, b

Scopigno, T., Ruocco, G., Sette, F., and Monaco, G.: Is the fragility of a liquid embedded in the properties of its glass?, Science, 302, 849–852, 2003. a

Segre, P. N., Liu, F., Umbanhowar, P., and Weitz, D. A.: An effective gravitational temperature for sedimentation, Nature, 409, 594–597, 2001. a

Seropian, G., Rust, A. C., and Sparks, R. S. J.: The gravitational stability of lenses in magma mushes: confined Rayleigh-Taylor instabilities, J. Geophys. Res.-Sol. Ea., 123, 3593–3607,, 2018. a

Shimomura, Y., Nishimura, T., and Sato, H.: Bubble growth processes in magma surrounded by an elastic medium, J. Volcanol. Geoth. Res., 155, 307–322,, 2006. a

Sigurdsson, H., Houghton, B., McNutt, S. R., Rymer, H., and Stix, J.: The Encyclopedia of Volcanoes, Elsevier,, 2015. a

Sparks, R.: The dynamics of bubble formation and growth in magmas: A review and analysis, J. Volcanol. Geoth. Res., 3, 1–37,, 1978. a

Sparks, R. S. J.: Forecasting volcanic eruptions, Earth Planet. Sc. Lett., 210, 1–15,, 2003. a

Stadtke, H.: Gasdynamic Aspects of Two-Phase Flow, 1st edn., Wiley-VCH, 288 pp., ISBN-10 352740578X, ISBN-13 978-3527405787, 2006. a, b, c, d, e, f, g, h

Suckale, J., Hager, B. H., Elkins-Tanton, L. T., and Nave, J.-C.: It takes three to tango: 2. Bubble dynamics in basaltic volcanoes and ramifications for modeling normal Strombolian activity, J. Geophys. Res.-Sol. Ea., 115, 2010a. a, b, c, d, e, f, g, h

Suckale, J., Nave, J.-C., and Hager, B. H.: It takes three to tango: 1. Simulating buoyancy-driven flow in the presence of large viscosity contrasts, J. Geophys. Res.-Sol. Ea., 115, B7,, 2010b. a

Suzuki, Y. J., Koyaguchi, T., Ogawa, M., and Hachisu, I.: A numerical study of turbulent mixing in eruption clouds using a three-dimensional fluid dynamics model, J. Geophys. Res.-Sol. Ea., 110, B8,, 2005. a, b

Thummala, P. P.: ReactingTwoPhaseEulerFoam description, Proceedings of CFD with OpenSource Software, (last access: 27 April 2022), 2016. a

Ubbink, O. and Issa, R.: A method for capturing sharp fluid interfaces on arbitrary meshes, J. Computat. Phys., 153, 26–50, 1999. a

Venier, C. M., Marquez Damian, S., and Nigro, N. M.: Numerical aspects of Eulerian gas–particles flow formulations, Comput. Fluids, 133, 151–169,, 2016. a

Wardle, K. E. and Weller, H. G.: Hybrid multiphase CFD solver for coupled dispersed/segregated flows in liquid-liquid extraction, Int. J. Chem. Eng., 2013, 128936,, 2013. a

Wark, D., Hildreth, W., Spear, F., Cherniak, D., and Watson, E.: Pre-eruption recharge of the Bishop magma system, Geology, 35, 235,, 2007. a

Winden, B.: An Open-Source Framework for Ship Performance CFD, in: SNAME 26th Offshore Symposium, OnePetro, 7 April 2021,, 2021. a

Xie, C., Tao, J., and Li, J.: Viscous Rayleigh-Taylor instability with and without diffusion effect, Appl. Math. Mech., 38, 263–270, 2017. a, b

Yeoh, G. H. and Tu, J. (Eds.): Computational techniques for multiphase flows, 2nd edn., Butterworth-Heinemann, 567–596,, 2019a. a, b, c

Yeoh, G. H. and Tu, J.: Chapter 10 - Future Trends in Handling Turbulent Multiphase Flows, in: Computational Techniques for Multiphase Flows (Second Edition), edited by: Yeoh, G. H. and Tu, J., Butterworth-Heinemann, second edn., 567–596,, 2019b. a

Short summary
Computer simulations play a fundamental role in understanding volcanic phenomena. The growing complexity of these simulations requires the development of flexible computational tools that can easily switch between sub-models and solution techniques as well as optimizations. MagmaFOAM is a newly developed library that allows for maximum flexibility for solving multiphase volcanic flows and promotes collaborative work for in-house and community model development, testing, and comparison.