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

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

Figure 1 Schematic representation of the (a) multiscale nature of a gas–liquid–solid flow (reprinted with minor modifications from Yeoh and Tu, 2019b, 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 Marschall, 2011).

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 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 2 MagmaFOAM–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 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₂O–CO₂ 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):

$$\begin{array}{}\text{(1)}& \mathit{\rho }(p,T,X)={\displaystyle \frac{M\left(X\right)}{V(p,T,X)}}={\displaystyle \frac{{\sum }_i{X}_i{M}_i}{{\sum }_i{X}_i{V}_i(p,T)}},\end{array}$$

where X_i is the mole fraction of the ith 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:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{V}_i(p,T)& =\sum _{l,m}{a}_{l,m}^i{T}^l{p}^m\\ & =a_{\mathrm{0},\mathrm{0}}^i+a_{\mathrm{1},\mathrm{0}}^{i}T+a_{\mathrm{0},\mathrm{1}}^{i}p+a_{\mathrm{1},\mathrm{1}}^{i}pT+\mathrm{\dots }\end{array}\end{array}$$

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₂O and CO₂ 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 η:

$$\begin{array}{}\text{(3)}& \log \mathit{\eta }=A+{\displaystyle \frac{B}{T-C}},\end{array}$$

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² s)⁻¹] and the interfacial area concentration A [m² m⁻³].

$$\begin{array}{}\text{(4)}& {\mathrm{\Gamma }}_i=J_{i}A\end{array}$$

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:

$$\begin{array}{}\text{(5)}& J_i=k_i\mathrm{\Delta }{C}_i,\end{array}$$

where k_i [kg (m² s)⁻¹] 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}}$):

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }{C}_i=C_i-C_i^{\mathrm{f}}.\end{array}$$

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

$$\begin{array}{}\text{(7)}& C_i^{\mathrm{f}}=C_i^{\mathrm{sat}}\left(p,T_{\mathrm{f}},X,X_{V\mathrm{tot}}\right),\end{array}$$

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₂O and CO₂, $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₂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.

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

Figure 3 Comparison 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 4 Rayleigh–Taylor instability ($A=\mathrm{0.5},\mathit{Re}=\mathrm{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 5 Rayleigh–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 (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).

Figure 6 Rayleigh–Taylor instability between a volatile-rich ($X_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}=\mathrm{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 (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}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ and $\mathit{\mu }={\mathrm{10}}^{-\mathrm{3}}\mathrm{Pa}\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⁻³. Surface tension is σ=0.073 N m⁻¹ for water, while a reasonable value for magmas is σ=0.15 N m⁻¹ (Colucci et al., 2016). In our tests we set σ=0.3 N m⁻¹ and ρ=3500 kg m⁻³ to be consistent with Suckale et al. (2010a). With a being the bubble radius and v₀ 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.

Figure 7 Simulations 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, Fr²≈0.4, We≈90, and $\mathrm{\Pi }={\mathrm{10}}^{-\mathrm{6}}$), (b) gradual breakup (Re≈25, Fr²≈0.3, We≈800, and $\mathrm{\Pi }={\mathrm{10}}^{-\mathrm{6}}$), and (c) catastrophic breakup (Re≈250, Fr²≈0.16, We≈1350, and $\mathrm{\Pi }={\mathrm{10}}^{-\mathrm{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 (${\mathit{Fr}}^{\mathit{2}}=\frac{{v_{\mathrm{0}}}^{\mathrm{2}}}{ga}$).

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₂ (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 8 Temporal evolution of bubble radius for an instantaneous decompression from p₀=150 MPa to p_L=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 H₂O. The red lines represent the same simulations with 1 wt % of CO₂ added in the melt. The diffusion coefficient of CO₂ is 1 order of magnitude smaller than H₂O. Initial conditions and parameters (see Appendix B) in all simulations are ρ_L=2300 kg m⁻³, $\mathit{\mu }=\mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}$ Pa s, σ=0.06 N m⁻¹, T=1123 K, $p_{\mathrm{G}}(t=\mathrm{0})=p_{\mathrm{0}}+\mathrm{2}\mathit{\sigma }/R(t=\mathrm{0})$, $R(t=\mathrm{0})={\mathrm{10}}^{-\mathrm{7}}$ m, $S_{\mathrm{0}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ m, and $C_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}^{\mathrm{0}}=\mathrm{0.053}\mathrm{wt}\mathrm{\%}$. 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 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.

Figure 9 Results at time t=0.015 s for the air–water shock tube using the SPWAT EOS. Solid lines: OpenFOAM; dashed lines: analytical solution (Stadtke, 2006); black lines: mixture; blue lines: liquid (water); red lines: gas (air). Mixture density is calculated a posteriori as ${\mathit{\rho }}_{\mathrm{mix}}={\mathit{\alpha }}_{\mathrm{g}}{\mathit{\rho }}_{\mathrm{g}}+(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{g}}){\mathit{\rho }}_{\mathrm{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: P_l=0.5 MPa, P_r=0.1 MPa, and $T_{\mathrm{l}}=T_{\mathrm{r}}=\mathrm{300}$ K for both phases; gas volume fraction: α_l=0.05, α_r=0.1, and $U_{\mathrm{l}}=U_{\mathrm{r}}=\mathrm{0}$ for both phases. Isobaric heat capacities of gas air and liquid water are Cp_g=1004.5 J kg⁻¹ K⁻¹ and Cp_l=4195 J kg⁻¹ K⁻¹, respectively (https://webbook.nist.gov/, 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⁻¹ m⁻¹ (https://webbook.nist.gov/, last access: 27 April 2022).

.

Figure 10 Results at time t=0.0065 s for the axisymmetric gas water–basalt shock tube using the Lange–Carmichael EOS (Lange and Carmichael, 1987) 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: P_l=10 MPa, P_r=0.1 MPa, and $T_{\mathrm{l}}=T_{\mathrm{r}}=\mathrm{1273}$ K for both phases; gas volume fraction: α_l=0.05, α_r=0.1, and $U_{\mathrm{l}}=U_{\mathrm{r}}=\mathrm{0}$ for both phases. Isobaric specific heat capacities in the gas and liquid phase are $C_{{\mathrm{P}}_{\mathrm{g}}}=\mathrm{2510}$ J kg⁻¹ K⁻¹ (https://webbook.nist.gov/, last access: 27 April 2022) as well as $C_{{\mathrm{P}}_{\mathrm{l},{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}=\mathrm{2278}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ and $C_{{\mathrm{P}}_{\mathrm{l},\mathrm{basalt}}}=\mathrm{1600}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, respectively (Lesher and Spera, 2015). Thermal conductivity of the liquid is 1.5 W K⁻¹ m⁻¹ (Lesher and Spera, 2015); for the water gas phase a Prandtl number of 0.9 is used, corresponding to a thermal conductivity of 0.14 W K⁻¹ m⁻¹ (https://webbook.nist.gov/, last access: 27 April 2022).

Figure 11 Results at time t=0.001 s for the axisymmetric gas water–basalt shock tube using the Lange–Carmichael EOS (Lange and Carmichael, 1987) 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: P_l=10 MPa, P_r=0.1 MPa, and $T_{\mathrm{l}}=T_{\mathrm{r}}=\mathrm{1273}$ K for both phases; gas volume fraction: α_l=0.05, α_r=1, and $U_{\mathrm{l}}=U_{\mathrm{r}}=\mathrm{0}$ for both phases. Isobaric specific heat capacities in the gas and liquid phase are $C_{{\mathrm{P}}_{\mathrm{g}}}=\mathrm{2510}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ (https://webbook.nist.gov/, last access: 27 April 2022) as well as $C_{{\mathrm{P}}_{\mathrm{l},{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}=\mathrm{2278}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ and $C_{{\mathrm{P}}_{\mathrm{l},\mathrm{basalt}}}=\mathrm{1600}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, respectively (Lesher and Spera, 2015). Thermal conductivity of the liquid is $\mathrm{1.5}\mathrm{W}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$ (Lesher and Spera, 2015); for the water gas phase a Prandtl number of 0.9 is used, corresponding to thermal conductivity of $\mathrm{0.14}W{K}^{-\mathrm{1}}{m}^{-\mathrm{1}}$ (https://webbook.nist.gov/, 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 10⁹ 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}}\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⁻¹ to 10⁻³ 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 12 Air–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: P_l=0.5 MPa, P_r=0.1 MPa, and $T_{\mathrm{l}}=T_{\mathrm{r}}=\mathrm{300}K$ for both phases; gas volume fraction: α_l=0.05, α_r=0.1, and $U_{\mathrm{l}}=U_{\mathrm{r}}=\mathrm{0}$ for both phases. Isobaric heat capacities of gas air and liquid water are $C_{{\mathrm{P}}_{\mathrm{g}}}=\mathrm{1004.5}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ and $C_{{\mathrm{P}}_{\mathrm{l}}}=\mathrm{4195}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, respectively (https://webbook.nist.gov/, 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⁻¹ m⁻¹, respectively (https://webbook.nist.gov/, 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., 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).

Figure 13 Reacting box simulation. At time 0 a small amount of gas (volume fraction ${\mathit{\alpha }}_{\mathrm{g}}={\mathrm{10}}^{-\mathrm{4}}$) is uniformly distributed in the box, and the basaltic melt is oversaturated in H₂O (5 wt %) and CO₂ (5000 ppm) at 80 MPa and 1373 K. The diffusion coefficients for H₂O and CO₂ in the basalt are 10⁻⁹ and ${\mathrm{10}}^{-\mathrm{10}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, respectively (Baker et al., 2005). Isobaric specific heat capacities in the gas phase are $C_{{\mathrm{P}}_{\mathrm{g},{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}=\mathrm{2900}Jk{g}^{-\mathrm{1}}{K}^{-\mathrm{1}}$ and $C_{{\mathrm{P}}_{\mathrm{g},{\mathrm{CO}}_{\mathrm{2}}}}=\mathrm{1390}Jk{g}^{-\mathrm{1}}{K}^{-\mathrm{1}}$ (https://webbook.nist.gov/, last access: 27 April 2022, Beaton et al., 1987) and in the liquid phase $C_{{\mathrm{P}}_{\mathrm{l},{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}=\mathrm{2278}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, $C_{{\mathrm{P}}_{\mathrm{l},{\mathrm{CO}}_{\mathrm{2}}}}=\mathrm{1600}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, and $C_{{\mathrm{P}}_{\mathrm{b}asalt}}=\mathrm{1600}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ (Lesher and Spera, 2015). Thermal conductivity of the liquid is $\mathrm{1.5}\mathrm{W}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$ (Lesher and Spera, 2015); for the gas phase Prandtl numbers of 0.9 for H₂O and 0.7 for CO₂ are used, corresponding to thermal conductivities of 0.16 and 0.09 $\mathrm{W}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$ (https://webbook.nist.gov/, 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 H₂O and CO₂. After $\approx \mathrm{4}\times {\mathrm{10}}^{\mathrm{5}}\mathrm{s}$, H₂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₂O (diffusion coefficient $D={\mathrm{10}}^{-\mathrm{9}}{\mathrm{m}}^{\mathrm{2}}{\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₂ is still out of thermodynamic equilibrium, as expected, because the diffusion coefficient is lower, being set to $D={\mathrm{10}}^{-\mathrm{10}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (Baker et al., 2005). The density and viscosity of the liquid vary with the decreasing dissolved H₂O. The gas density decreases because of increasing H₂O and decreasing CO₂ 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 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).

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 (Co_max=0.01).

Figure A1 Time evolution of the amplitude of two single-mode perturbations (k=0.5 – a; k=0.15 – b) 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 A2 Extrapolated 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.

$$\begin{array}{ll}\text{(B1)}& {\displaystyle }& {\displaystyle }{\mathit{\rho }}_{\mathrm{L}}R{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}R}{\mathrm{d}{t}^{\mathrm{2}}}}\left(\mathrm{1}-{\displaystyle \frac{R}{S}}\right)+{\mathit{\rho }}_{\mathrm{L}}{\left({\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}}\right)}^{\mathrm{2}}\left({\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}-{\displaystyle \frac{\mathrm{2}R}{S}}+{\displaystyle \frac{R^{\mathrm{4}}}{\mathrm{2}{S}^{\mathrm{4}}}}\right)={\displaystyle }& {\displaystyle }{p}_{\mathrm{G}}-p_{\mathrm{L}}+\mathrm{4}{\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}}{R}^{\mathrm{2}}\left(-\mathrm{3}\underset{R}{\overset{S}{\int }}{\displaystyle \frac{\mathit{\mu }\left(r\right)}{r^{\mathrm{4}}}}\mathrm{d}r\right)-{\displaystyle \frac{\mathrm{2}\mathit{\sigma }}{R}}\end{array}$$

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

$$\begin{array}{}\text{(B2)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left(R^{\mathrm{3}}{\mathit{\rho }}_{\mathrm{G}}\right)=\mathrm{3}{R}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{L}}\sum _{i=\mathrm{1}}^N{D}_i{\left[{\displaystyle \frac{\partial C_i}{\partial r}}\right]}_{r=R},\end{array}$$

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

$$\begin{array}{}\text{(B3)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left(R^{\mathrm{3}}{\mathit{\rho }}_{\mathrm{G}}{Y}_i\right)=\mathrm{3}{R}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{L}}{D}_i{\left[{\displaystyle \frac{\partial C_i}{\partial r}}\right]}_{r=R},\end{array}$$

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

$$\begin{array}{}\text{(B4)}& \begin{array}{r}{\left[{\displaystyle \frac{\mathrm{d}{C}_i}{\mathrm{d}r}}\right]}_{r=R}\\ ={\displaystyle \frac{S_{\mathrm{0}}^{\mathrm{3}}\left(C_i^{\mathrm{0}}-C_i^{\mathrm{sat}}\right)+\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{L}}}\left({\left[R^{\mathrm{3}}{\mathit{\rho }}_{\mathrm{G}}{Y}_i\right]}_{t=\mathrm{0}}-R^{\mathrm{3}}{\mathit{\rho }}_{\mathrm{G}}{Y}_i\right)}{S_{\mathrm{0}}^{\mathrm{3}}R-\frac{\mathrm{3}}{\mathrm{2}}\left(S^{\mathrm{2}}-R^{\mathrm{2}}\right)R^{\mathrm{2}}}},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B5)}& \left(-\mathrm{3}\underset{R}{\overset{S}{\int }}{\displaystyle \frac{\mathit{\mu }\left(C_i\right(r\left)\right)}{r^{\mathrm{4}}}}\mathrm{d}r\right)=\mathit{\mu }\left({\displaystyle \frac{\mathrm{1}}{S^{\mathrm{3}}}}-{\displaystyle \frac{\mathrm{1}}{R^{\mathrm{3}}}}\right).\end{array}$$

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

$$\begin{array}{}\text{(B6)}& \mathit{\alpha }={\displaystyle \frac{R^{\mathrm{3}}}{S^{\mathrm{3}}}}.\end{array}$$

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 C1 Results 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: P_l=0.1 MPa, P_r=0.01 MPa; T_l=348.432 K, T_r=278.746 K; gas volume fraction α_l=1, α_r=0; $U_{\mathrm{l}}=U_{\mathrm{r}}=\mathrm{0}$. Isobaric heat capacity is $C_{\mathrm{P}}=\mathrm{1004.5}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, corresponding to heat capacity ratio γ=1.4.

Figure C2 Results 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: P_l=10 MPa, P_r=0.1 MPa; $T_{\mathrm{l}}=T_{\mathrm{r}}=\mathrm{300}\mathrm{K}$; gas volume fraction ${\mathit{\alpha }}_{\mathrm{l}}={\mathit{\alpha }}_{\mathrm{r}}=\mathrm{0}$; $U_{\mathrm{l}}=U_{\mathrm{r}}=\mathrm{0}$. Isobaric heat capacity is $C_{\mathrm{P}}=\mathrm{4195}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ (https://webbook.nist.gov/, last access: 27 April 2022).

Figure C3 Results 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: P_l=10 MPa, P_r=0.1 MPa; $T_{\mathrm{l}}=T_{\mathrm{r}}=\mathrm{1373}\mathrm{K}$; gas volume fraction ${\mathit{\alpha }}_{\mathrm{l}}={\mathit{\alpha }}_{\mathrm{r}}=\mathrm{0}$; $U_{\mathrm{l}}=U_{\mathrm{r}}=\mathrm{0}$. Isobaric heat capacity is $C_{\mathrm{P}}=\mathrm{1600}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ (Lesher and Spera, 2015). Thermal conductivity is $\mathrm{1.5}\mathrm{W}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$ (Lesher and Spera, 2015).

Figure C4 Results 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: P_l=0.5 MPa, P_r=0.1 MPa; $T_{\mathrm{l}}=T_{\mathrm{r}}=\mathrm{300}\mathrm{K}$ for both phases; gas volume fraction α_l=0.05, α_r=1; $U_{\mathrm{l}}=U_{\mathrm{r}}=\mathrm{0}$ for both phases. Isobaric heat capacities of gas air and liquid water are $C_{{\mathrm{P}}_{\mathrm{g}}}=\mathrm{1004.5}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ and $C_{{\mathrm{P}}_{\mathrm{l}}}=\mathrm{4195}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, respectively (https://webbook.nist.gov/, 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 $\mathrm{W}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$ (https://webbook.nist.gov/, 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 D1 Oxide compositions for the magmas used in benchmarking simulations. Amounts are relative.