**Model description paper**| 10 May 2022

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

Federico Brogi Simone Colucci Jacopo Matrone Chiara Paola Montagna Mattia De' Michieli Vitturi and Paolo Papale

^{1,2},

^{1},

^{3},

^{1},

^{1,4},

^{1}

**Federico Brogi et al.**Federico Brogi Simone Colucci Jacopo Matrone Chiara Paola Montagna Mattia De' Michieli Vitturi and Paolo Papale

^{1,2},

^{1},

^{3},

^{1},

^{1,4},

^{1}

^{1}Istituto Nazionale di Geofisica e Vulcanologia, sezione di Pisa, Pisa, Italy^{2}Istituto Nazionale di Oceanografia e di Geofisica Sperimentale, Trieste, Italy^{3}Dipartimento di Matematica, Università degli Studi di Firenze, Florence, Italy^{4}Department of Geology, University at Buffalo, Buffalo, New York, USA

^{1}Istituto Nazionale di Geofisica e Vulcanologia, sezione di Pisa, Pisa, Italy^{2}Istituto Nazionale di Oceanografia e di Geofisica Sperimentale, Trieste, Italy^{3}Dipartimento di Matematica, Università degli Studi di Firenze, Florence, Italy^{4}Department of Geology, University at Buffalo, Buffalo, New York, USA

**Correspondence**: Federico Brogi (federico.brogi@ingv.it)

**Correspondence**: Federico Brogi (federico.brogi@ingv.it)

Received: 28 Jun 2021 – Discussion started: 12 Aug 2021 – Revised: 16 Mar 2022 – Accepted: 24 Mar 2022 – Published: 10 May 2022

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.

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 (Sparks, 2003; 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 (Melnik, 2000; Papale, 2001; 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; Dufek, 2016), and lava flows (Griffiths, 2000). 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 Koyaguchi, 2005; Macedonio et al., 2005; Sahagian and Proussevitch, 2005; 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 Suckale, 2019).

Volcanic transport processes are typically characterized by a wide range of spatial and temporal scales at which different interacting physical subprocesses occur (Griffiths, 2000; Gonnermann and Manga, 2007; Dufek, 2016). 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.

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 Mahesh, 1998), 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 Suckale, 2020) 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 Suckale, 2020). 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 Tu, 2019a). 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 Suckale, 2019; 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 (Marschall, 2011). 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 Fox, 2007), including magmas (Keller and Suckale, 2019). 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 (Papale, 2001; Dufek and Bergantz, 2005), 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 Balachandar, 2001). 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 Sparks, 2006) 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 Brennen, 2005). 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 Tu, 2019a).

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., Winden, 2021). Moreover, the recently established exaFoam consortium will improve computational performance, enabling the “OpenFOAM community” to efficiently exploit “the current evolving HPC hardware and middleware” (http://www.exafoam.eu, 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 Carmichael, 1987) 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.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.

### 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 *p*–*T* 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 Spera, 2015).
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 H_{2}O–CO_{2} saturation over a broad range of *p*–*T* 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 (*M*_{i}) and molar volumes (*V*_{i}):

where *X*_{i} is the mole fraction of the *i*th oxide component. To a good approximation, molar volumes do not depend on melt composition (Lesher and Spera, 2015) and can be computed with a polynomial expansion:

The polynomial coefficients ${a}_{l,m}^{i}$ 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 H_{2}O and CO_{2} 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 *η*:

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
($\sim {\mathrm{10}}^{-\mathrm{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 *T*≤*C* (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 *J*_{i} [kg (m^{2} s)^{−1}] and the interfacial area concentration *A* [m^{2} m^{−3}].

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 Hibiki, 2006).
The model for *J*_{i} expresses the driving force for diffusive mass transfer of the component *i* and can be calculated with the following relationship:

where *k*_{i} [kg (m^{2} s)^{−1}] is the mass transfer coefficient, a function of the diffusion coefficient *D*_{i} (Cussler, 2009; Thummala, 2016), and
Δ*C*_{i} is the difference
between the mass fraction of the species in the phase (*C*_{i}) and at the interface (${C}_{i}^{\mathrm{f}}$):

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

where ${C}_{i}^{\mathrm{sat}}$ 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,(*T*_{f}), melt composition (*X*), and total amount of
volatiles of all species (*X*_{Vtot}).
For magmas with H_{2}O and CO_{2}, ${C}_{i}^{\mathrm{sat}}$ can be computed using SOLWCAD (Papale et al., 2006) or other dedicated models (e.g., Newman and Lowenstern, 2002; 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 ${C}_{i}^{\mathrm{sat}}$ 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 Sahagian, 1998). 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 Sahagian, 1998; 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 (H_{2}O) 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 (Manga, 1996).
All model equations can be found in Appendix B and are solved as systems of ordinary differential equations (ODEs) using the OpenFOAM ODE solvers.

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 Kothe, 1998; 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 Issa, 1999).
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 Poli, 2012; 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 Bergantz, 2003; 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
(Chandrasekhar, 1955), surface tension
(Chandrasekhar, 2013), compressibility (Mitchner and Landshoff, 1964), and
diffusivity (Duff et al., 1962; Xie et al., 2017). The problem parameters can be
expressed by two dimensionless numbers: the Atwood number
$\mathrm{Atw}=({\mathit{\rho}}_{h}-{\mathit{\rho}}_{l})/({\mathit{\rho}}_{h}+{\mathit{\rho}}_{l})$ and the Reynolds
number ($\mathit{Re}=\sqrt{Wg}W/\mathit{\nu}$), 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\times \mathrm{2}W$). 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=\mathrm{2}\mathit{\pi}/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.

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 (*Co*_{max}=0.5) to speed up the simulation. For a given mesh resolution the accuracy and convergence of the solution depend on the values of *Co*_{max} and the number of iterations (`nIter`

) used to solve the pressure–velocity coupling with the PISO algorithm (Issa, 1986). Generally, larger *Co*_{max} (<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 *Co*_{max}, 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 (${X}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}=\mathrm{2}$ wt %) and a chemically more
evolved andesitic melt. Temperature
is set to *T*=1300 ^{∘}C 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 *Co*_{max}=0.5). As a result, the simulation is computationally much more demanding. The simulated time covers the entire overturning process (Fig. 6).

### 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
(Grace, 1973). 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 $\mathit{\rho}={\mathrm{10}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{3}}$
and $\mathit{\mu}={\mathrm{10}}^{-\mathrm{3}}\phantom{\rule{0.125em}{0ex}}\mathrm{Pa}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$, 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 *v*_{0} 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 $\mathit{Re}=\frac{\mathit{\rho}{v}_{\mathrm{0}}a}{\mathit{\mu}}$ (inertial to viscous forces), Froude number $\mathit{Fr}=\frac{{v}_{\mathrm{0}}}{\sqrt{ga}}$ (inertia to buoyancy forces), Weber number $\mathit{We}=\frac{\mathit{\rho}{v}_{\mathrm{0}}^{\mathrm{2}}a}{\mathit{\sigma}}$ (inertia to surface tension), and gas to liquid viscosity ratio $\mathrm{\Pi}={\mathit{\mu}}_{g}/\mathit{\mu}$. 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 $\mathrm{\Pi}={\mathrm{10}}^{-\mathrm{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.

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 CO_{2} (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).

## 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 Manga, 2007). 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 Mitani, 2005; 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 (Stadtke, 2006). 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 (Stadtke, 2006). 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 (Stadtke, 2006). 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 (Stadtke, 2006)
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 (Stadtke, 2006). 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 (Stadtke, 2006) (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 (Stadtke, 2006) 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 ($\mathrm{0.5}<\mathit{\alpha}<\mathrm{0.7}$; e.g., Sparks, 1978, 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.

*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 10^{9} Pa 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 (${\mathit{\tau}}_{\mathrm{b}}\propto {D}_{\mathrm{b}}^{\mathrm{2}}/{\mathit{\nu}}_{\mathrm{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 (${\mathit{\tau}}_{\mathrm{b}}\approx {\mathrm{10}}^{-\mathrm{6}}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$), 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 Lo, 1999). The Partial Elimination Algorithm (PEA), implemented in OpenFOAM, shows the best convergence performance compared to other algorithms (Karema and Lo, 1999; 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.

### 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., *k*_{i} 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 Hibiki, 2006) 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 Hibiki, 2006).

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 H_{2}O and CO_{2}. After
$\approx \mathrm{4}\times {\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$, H_{2}O 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 H_{2}O
(diffusion coefficient $D={\mathrm{10}}^{-\mathrm{9}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, Baker et al., 2005) around a bubble
with radius *R*≈2 cm (${\mathit{\tau}}_{\mathrm{d}}={R}^{\mathrm{2}}/D$;
Lensky et al., 2004). The dissolved
CO_{2} is still out of thermodynamic equilibrium, as expected,
because the diffusion coefficient is lower, being set to $D={\mathrm{10}}^{-\mathrm{10}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ (Baker et al., 2005). The density and viscosity of the liquid vary with the decreasing dissolved H_{2}O. The gas density decreases because of increasing H_{2}O and decreasing CO_{2} that produce a decrease in the molar mass of the gas mixture.

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 Fox, 2013) will improve our understanding of how polydispersity can impact magmatic system evolution (Colucci et al., 2017a; de' Michieli Vitturi and Pardini, 2021). 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 Tu, 2019a). 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 Fox, 2013). 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 Weller, 2013). From a volcanological perspective, predicting flow regime changes is of crucial importance to understand effusive–explosive transitions in eruptive activity (Gonnermann and Manga, 2007).

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 (*Co*_{max}=0.01).

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 (${S}^{\mathrm{3}}={S}_{\mathrm{0}}^{\mathrm{3}}+{R}^{\mathrm{3}}$; ${S}_{\mathrm{0}}^{\mathrm{3}}={S}^{\mathrm{3}}(t=\mathrm{0})-{R}^{\mathrm{3}}(t=\mathrm{0})$), *p*_{G} is the gas pressure inside the bubble, *p*_{L} 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 *p*_{L}(*t*) this represents an equation that can be solved to find *R*(*t*) provided *p*_{G}(*t*). *p*_{G} 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

where *D*_{i} is the mass diffusivity and *C*_{i} the concentration of the *i*th species dissolved in the melt.
Mass conservation of the *i*th dissolved species is given by

where *Y*_{i} is the concentration in the gas phase.
Assuming local thermodynamic equilibrium at the bubble–melt interface
(i.e., ${C}_{i}\left(R\right)={C}_{i}^{\mathrm{sat}}\left({p}_{\mathrm{G}}\right)$, where ${C}_{i}^{\mathrm{sat}}$ 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

where ${C}_{i}^{\mathrm{0}}$ is the concentration in the melt at time 0. Assuming constant viscosity, the term in Eq. (B1) is analytically integrated to obtain

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

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

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 (https://doi.org/10.5281/zenodo.5031825, Brogi et al., 2021).

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.

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.

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).

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, https://doi.org/10.1093/petrology/egi084, 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, https://doi.org/10.1002/2016GL071875, 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, https://doi.org/10.1038/ngeo2534, 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], https://doi.org/10.5281/zenodo.5031825, 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, https://doi.org/10.5194/gmd-9-697-2016, 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, https://doi.org/10.1017/S0305004100030048, 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, https://doi.org/10.1029/2005JB004174, 2006. a

Colucci, S., Battaglia, M., and Trigila, R.: A thermodynamical model for the
surface tension of silicate melts in contact with *H*_{2}*O* gas, Geochim.
Cosmochim. Ac., 175, 127, https://doi.org/10.1016/j.gca.2015.10.037, 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, https://doi.org/10.1007/s00410-017-1421-6, 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, https://doi.org/10.1002/2016JB013383, 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, https://doi.org/10.1016/j.jvolgeores.2020.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, https://doi.org/10.1017/CBO9780511805134, 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, https://doi.org/10.5194/gmd-14-1345-2021, 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, https://doi.org/10.1016/j.epsl.2008.05.025, 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, https://doi.org/10.5194/gmd-8-2447-2015, 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, https://doi.org/10.1088/1749-4699/5/1/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, https://doi.org/10.1038/nature10706, 2012. a

Dufek, J.: The Fluid Mechanics of Pyroclastic Density Currents, Annu. Rev. Fluid Mech., 48, 459–485, https://doi.org/10.1146/annurev-fluid-122414-034252, 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, https://doi.org/10.1016/j.parco.2007.04.003, 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, https://doi.org/10.1038/s41598-019-40160-1, 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, https://doi.org/10.1016/j.ijmultiphaseflow.2018.10.010, 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, https://doi.org/10.1146/annurev.fluid.39.050905.110207, 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, https://doi.org/10.1146/annurev.fluid.32.1.477, 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, https://doi.org/10.1038/s41586-018-0746-2, 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, https://doi.org/10.1029/1998JB900116, 1999. a

Karema, H. and Lo, S.: Efficiency of interphase coupling algorithms in fluidized bed conditions, Comput. Fluids, 28, 323–360, https://doi.org/10.1016/S0045-7930(98)00028-0, 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, https://doi.org/10.1093/gji/ggz287, 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, https://doi.org/10.1029/2004JB003513, 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, https://doi.org/10.1016/B978-0-12-385938-9.00005-5, 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, https://doi.org/10.1007/s004450050122, 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, https://doi.org/10.1016/j.jvolgeores.2004.09.015, 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, https://doi.org/10.1007/978-3-211-72464-4, 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, https://doi.org/10.1017/CBO9781139016599, 2013. a, b

Marschall, H.: Towards the numerical simulation of multi-scale two-phase flows, PhD thesis, Technische Universität München, https://nbn-resolving.org/urn/resolver.pl?urn:nbn:de:bvb:91-diss-20111222-1080878-1-7, 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, https://doi.org/10.1016/j.jvolgeores.2020.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, https://doi.org/10.1016/j.jvolgeores.2004.09.011, 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, https://doi.org/10.1073/pnas.0911705106, 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, https://doi.org/10.1007/s004450000072, 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, https://doi.org/10.1017/jfm.2019.276, 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, https://doi.org/10.1007/11157_2017_30, 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–H_{2}O–CO_{2}
solution model written in Visual Basic for excel, Comput. Geosci.,
28, 597–604, https://doi.org/10.1016/S0098-3004(01)00081-4, 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, https://doi.org/10.1029/2000JB900428, 2001. a, b

Papale, P., Moretti, R., and D., B.: The compositional dependence of the
saturation surface of H_{2}O+CO_{2} 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, https://doi.org/10.1016/j.lithos.2012.02.002, 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, https://doi.org/10.1007/s00445-009-0329-z, 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, https://doi.org/10.1029/2019JB018549, 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, https://doi.org/10.1016/j.ces.2011.02.030, 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, https://doi.org/10.1029/2008GC002022, 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, https://doi.org/10.1016/j.jvolgeores.2004.09.016, 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, https://www.researchgate.net/publication/280099401_Direct_Numerical_Simulation_of_Single_Bubble_Rising_in_Viscous_Stagnant_Liquid (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, https://doi.org/10.1029/2018JB015523, 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, https://doi.org/10.1016/j.jvolgeores.2006.04.003, 2006. a

Sigurdsson, H., Houghton, B., McNutt, S. R., Rymer, H., and Stix, J.: The Encyclopedia of Volcanoes, Elsevier, https://doi.org/10.1016/C2015-0-00175-7, 2015. a

Sparks, R.: The dynamics of bubble formation and growth in magmas: A review and analysis, J. Volcanol. Geoth. Res., 3, 1–37, https://doi.org/10.1016/0377-0273(78)90002-1, 1978. a

Sparks, R. S. J.: Forecasting volcanic eruptions, Earth Planet. Sc. Lett., 210, 1–15, https://doi.org/10.1016/S0012-821X(03)00124-9, 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, https://doi.org/10.1029/2009JB006916, 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, https://doi.org/10.1029/2004JB003460, 2005. a, b

Thummala, P. P.: ReactingTwoPhaseEulerFoam description, Proceedings of CFD with OpenSource Software, http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2016/PhanindraPrasadThummala/reactingTwoPhaseEulerFoam.pdf (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, https://doi.org/10.1016/j.compfluid.2016.05.003, 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, https://doi.org/10.1155/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, https://doi.org/10.1130/G23316A.1, 2007. a

Winden, B.: An Open-Source Framework for Ship Performance CFD, in: SNAME 26th Offshore Symposium, OnePetro, 7 April 2021, https://doi.org/10.5957/TOS-2021-22, 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, https://doi.org/10.1016/B978-0-08-102453-9.00010-6, 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, https://doi.org/10.1016/B978-0-08-102453-9.00010-6, 2019b. a

- Abstract
- Introduction
- MagmaFOAM ingredients
- Benchmarks and test cases
- Conclusions
- Appendix A: Linear Rayleigh–Taylor instability
- Appendix B: multiComponentODERPShellDStatic: model equations
- Appendix C: Shock Tube
- Appendix D: Magmatic compositions
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- MagmaFOAM ingredients
- Benchmarks and test cases
- Conclusions
- Appendix A: Linear Rayleigh–Taylor instability
- Appendix B: multiComponentODERPShellDStatic: model equations
- Appendix C: Shock Tube
- Appendix D: Magmatic compositions
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References