Understanding magma transport in the crust is one of the major challenges in modern volcanology. Current models depict magmatic systems as an interconnected network of compositionally heterogeneous magmas, involving multiple chambers and dikes that extend from shallow to deep levels in the crust (Gudmundsson, 1995; Gudmundsson, 2011; Blundy and Annen, 2016).

Magma chambers are the main engines of active volcanoes and serve as a storage region for magma ascent and its chemical evolution (DePaolo, 1981; Bachmann and Bergantz, 2003). The pressure in the chamber may vary largely over time due to a variety of processes, including fractional crystallization, volatile exsolution and magma recharge, potentially leading towards an eruption (Sparks and Huppert, 1984; Folch and Marti, 1998; Gudmundsson, 2012; Edmonds and Wallace, 2017; Papale et al., 2017; Cassidy et al., 2018). Therefore, studying magma dynamics in a chamber is of utmost importance and has become a mainstream theme over the recent years (Bergantz, 2000; Couch et al., 2001; Gutiérrez and Miguel, 2010; Karlstrom et al., 2010; Bergantz et al., 2015; Garg et al., 2019) as the ongoing flow dynamics and associated surface signals can depict the current state of the volcano and its possible evolutions (Gottsmann et al., 2011; Longo et al., 2012a; Bagagli et al., 2017; Sparks and Cashman, 2017; Carrara, 2019; Lieto et al., 2020).

Typically, buoyancy instabilities develop when a low-density, volatile-rich, primitive, hot magma ascends in the crust and interacts with previously stored, more chemically evolved, partially degassed and denser magma in shallow chambers (Semenov and Polyansky, 2017; Garg et al., 2019). There are multiple pieces of evidence of eruptions shortly preceded by interaction of compositionally different magmas at shallow depths (Tomiya et al., 2013; Sides et al., 2014; Perugini et al., 2015; Sundermeyer et al., 2020). The degree of magma mixing, which spans a continuum from mechanical mingling to complete chemical homogenisation, depends upon the magma properties, the driving forces, and the time available for mixing (Garg et al., 2019).

Recently a number of authors [e.g., (Sparks and Cashman, 2017), and references therein] have been hypothesizing that although usually not seen among the erupted volcanic products, a crystal-rich mush may constitute a large dominant portion of the magmatic plumbing system, with dominantly liquid magma possibly being an ephemeral occurrence driven by melt segregation and contributing to characterize the unrest and eruption states of the volcano. Clear evidence from active magmatic systems is not available yet, as volcanic plumbing systems continue to be hidden from direct observation. A few accidental encounters with buried shallow magma bodies during geothermal drilling do not appear to support the presence of important layers of mushy magma across the rock-magma transition (Eichelberger, 2020). More is needed before we get to a clear, robust picture of active magmatic systems, justifying substantial efforts to overcome the current limitations in directly observing, measuring and sampling underground magma (https://www.kmt.is). Here we focus on either dominantly liquid reservoirs, or on the dominantly liquid core region of mushy magma reservoirs, over time scales that are typical of individual magma convection events, likely much less than the time scales associated with the existence of ephemeral liquid-dominated magmatic bodies.

Over the years significant efforts have been invested to study several aspects of the magma physics relevant to understand the volcano dynamics and anticipate their evolution. The physics of magma mixing has been studied through numerical simulations and experiments with both synthetic and natural compositions (Campbell and Turner, 1986; Huppert et al., 1986; Jellinek and Kerr, 1999; Michioka and Sumita, 2005; Sato and Sato, 2009; De Campos et al., 2011; Laumonier et al., 2014; Morgavi et al., 2015; Perugini et al., 2015; Bergantz et al., 2015; Schleicher et al., 2016; Garg et al., 2019). Usually, numerical simulations are simplified by referring to a 2D geometrical configuration, because 3D dynamics can be very complex and are extremely expensive in terms of required computational resources. Even 2D magma mixing simulations need extremely refined meshes to capture the flow features down to the dm (or lower) scales that are sometimes required, e.g., at the intersection between feeding dykes and chambers (Longo et al., 2012a; Papale et al., 2017). Solving such problems in 3D on a single computer is practically impossible. Nowadays the simulations are usually run over clusters of cores providing high speed data processing capability. In high performance computing (HPC) paradigm, many processors work simultaneously to produce exceptional computational power and to significantly reduce the total computational time (Dowd and Charles, 2010; Flanagan et al., 2020). High performance is achieved through parallel computing in which large computational domains are subdivided into smaller, interconnected ones, which can be solved at the same time (Gottlieb, 1989; Asanovic et al., 2006).

So far the numerical simulations performed for magma mixing in the literature are in 2D (Bergantz, 2000; Longo et al., 2012a; Papale et al., 2017; Garg et al., 2019). However, strictly speaking, 2D calculations apply only to flows that are inherently two-dimensional, and can be extended to real 3D systems only under restricted conditions whereby any change in physical quantities over the third dimension can be neglected. Real magmatic systems are commonly such that 2D simplifications can only be introduced with some arbitrariness, typically with the objective of extracting zero-order approximations of more complex 3D processes and dynamics. In most cases, however, the loss of information due to 2D simplification is unknown.

The purpose of this work is that of 1) presenting a physical and numerical model for transient 3D magma dynamics, including its parallel implementation and scaling performance, 2) applying the model to 3D magma chamber convection dynamics in an initially stratified, gravitationally unstable magma chamber, and 3) comparing the 3D simulation results with the 2D case.

Natural magmas are composed of crystals, melt oxides and volatiles, with the latter that can be dissolved in the liquid phase or exsolved in a gas phase with proportions depending on local pressure (p), temperature (T), and composition (Y). Flow conditions span wide ranges from essentially incompressible to largely compressible, including supersonic flow during eruptions. In magmas where the volatiles are largely dissolved, and the Mach number is ≪ 1, the flow remains weakly compressible. In this work we are specifically interested in modelling the flow dynamics inside magma chambers where two compositionally different magmas interact with each other. The model that we present is however general, and can be applied to transonic flow as well.

We model compressible/incompressible flow of multi-component magma with GALES (Longo et al., 2012b; Garg et al., 2018a; Garg et al., 2018b; Garg et al., 2021). The numerical scheme and parallel implementation in GALES are described in Appendixes A, B. GALES has been validated on a number of 2D/3D benchmarks for multi-component compressible/incompressible flows (Longo et al., 2012b), single-component compressible/incompressible flows (Garg et al., 2018a), free surface flows (Garg et al., 2018b) and fluid-structure interaction (Garg et al., 2021). Additional validation tests for compressible and incompressible flows with 3D geometries are reported in the Supplementary Material.

The properties of each magma component in GALES are computed from local conditions including composition in terms of ten major oxides plus the two major volatile species H₂O and CO₂. The mixture is assumed to be in chemical, mechanical and thermodynamic equilibrium. The flow is governed by the following equations: ∂ ρ ∂ t + ∇ ⋅ ρ v = 0 (1) ∂ ρ v ∂ t + ∇ ⋅ ρ v ⊗ v = ∇ ⋅ σ + ρ b , (2) ∂ ρ E ∂ t + ∇ ⋅ ρ E v = ∇ ⋅ σ v − q − h ′ J ′ + ρ b ⋅ v , (3) ∂ ρ Y k ∂ t + ∇ ⋅ ρ v Y k = − ∇ ⋅ J k k = 1,2 , … , n (4)

Equations 1–3 represent mass, momentum and energy conservation of the mixture, while Eq. 4 is mass conservation of the kth component in the n-component mixture. The equations are same as in (Garg et al., 2019) with a three-dimensional extension of vector and tensor quantities. In the above equations, Y _k is the mass fraction of the kth component, with ∑ _k Y _k = 1. This implies that for given n components, only n − 1 components are independent and require each one expression of Eq. 4. For an n-component mixture, the density (ρ) is given by: 1 ρ = Y 1 ρ 1 + Y 2 ρ 2 + ⋯ + Y n ρ n (5) where ρ _k is the density of the kth component which depends on local p-T and composition; v is the velocity; b is body force vector per unit mass; E = c _v T + |v| ²/2 is total specific energy, with c _v being the specific heat capacity at constant volume; T is temperature; h is the vector of specific enthalpies for the components.

The mass diffusion flux is modelled with Fick’s law as J k = − ρ D ∇ Y k (6) where D is the mass diffusion coefficient matrix. Viscous flux is modelled as σ = μ ∇ v + ∇ v ′ − 2 3 μ ∇ ⋅ v I − p I = τ − p I (7) where p is pressure, μ is the viscosity of the mixture and τ is the viscous stress tensor. The heat flux is modelled by Fourier’s law: q = − κ ∇ T (8) where κ is thermal conductivity.

Equations 1–4 are written in the conservative form and describe fully compressible flow. The incompressible flow equations are merely a simplification of the above equations by referring to a constant density. We remark that although in magma reservoirs the Mach number is usually very low, the density of the magmatic mixture varies, mostly as a response to phase changes of volatile components. Therefore, considering the flow as fully incompressible would miss many important flow features of gas-bearing magmas.

In magma reservoirs, and over the time scale of individual convection events analysed here, phase separation is of minimum importance. In our previous works (Longo et al., 2012a; Papale et al., 2017; Garg et al., 2019) we have estimated that flow Stokes number St = t ₀ v ₀/l ₀, where t ₀ is the relaxation time of the crystal, v ₀ is the flow velocity and l ₀ is the diameter of the crystal, remains very low, of order 10^–4 for crystals up to cm size. The same is true for gas bubbles. Based on the typical gas volume fractions obtained in our previous works and assumed bubble number density as low as 10¹⁴ m⁻³ (Cashman, 2000), the value of St remains less than 10^–4. The low value of St indicates that mechanical phase separation is negligible, and the relative velocity terms describing phase separation can be safely ignored.

The physical properties of magma are modelled as a function of local pressure, temperature and composition in the space-time domain, as they evolve during magma convection and mixing. Throughout this paper the word “mixing” refers to the scale of our analysis, whereby the smallest elements in the computational domain have linear dimension of order 1 m. At such a scale many orders of magnitude larger than the scale of molecules, and for the employed computational times of order a few hours, chemical mixing is unresolved and likely to be poorly effective. Therefore, we refer only to mechanical mixing whereby both magma types are present within individual computational elements, with no reference to physical processes occurring at the molecular scale.

As components for use in Eq. 4 above, we refer here to the magma types involved in convection and mixing (e.g., “andesite”, “rhyolite”, etc.), each expressed in terms of oxides including volatile species. Accordingly, each component is modeled as a mixture of 1) silicate melt including the dissolved volatiles and 2) exsolved volatiles. Non-reactive crystals can be added and their effects in modifying the mechanical and thermal properties of magma can be accounted for (for simplicity, however, we have neglected crystals in the computations illustrated below). The density of the volatile free silicate melt is modelled according to (Lange, 1994) and the effects on the density by dissolved H₂O is computed by the model of (Burnham et al., 1969). For the gas phase, we use the ideal gas equation as the equation of state. For each magma, the phase distribution of volatiles is computed by multi-component (H₂O + CO₂) gas-melt equilibrium modelling (Papale et al., 2006). The viscosity of each magma component is modelled as a function of oxide composition, dissolved H₂O and temperature (non-Arrhenian) as described in (Giordano et al., 2008), with the effect of non-deformable gas bubbles accounted for by the model of (Ishii and Zuber, 1979) and strain-rate dependent non-Newtonian rheology due to the presence of crystals and of deforming bubbles (Caricchi et al., 2007; Pal, 2003). The viscosity of the two-component mixture is modelled as μ = ∑ _k x _k μ _k , where x _k and μ _k are, respectively, the mole fraction and the viscosity of the kth mixture component (again, in our case, “andesite”, “rhyolite”, etc.). The specific heat capacity at constant pressure for the melt is computed as c _p = ∑ _j x _j c _pj , where c _pj and x _j are, respectively, the specific heat capacity at constant pressure and mole fraction for the jth oxide subcomponent [Table 1 in (Garg et al., 2019)]. The specific heat capacity at constant volume, c _v is computed as c _v = c _p − α ² T/(βρ), where α and β are isobaric expansion coefficient and isothermal compressibility coefficient, respectively. In this work we use a constant thermal conductivity, κ = 1.2 Wm ⁻¹ K ⁻¹ (Garg et al., 2019).

The numerical scheme and the stabilization terms employed for the solution of the set of Eqs 1–4 above, are reported in the Appendix.

As discussed above, magmatic systems often consist of multiple heterogeneous magmas stored in a network of interconnected sills and dykes. Magma mixing is understood as one of many possible mechanisms for triggering an eruption (Sparks et al., 1977). Petrology and geochemistry analysis indicates that multiple intrusions of mafic to intermediate magma in more chemically evolved magmas stored at shallow depths may produce hybrid melts with zoned crystals. Sometimes the ejecta contain abundant inter-mingled hybrid magmas, suggesting efficient stirring and mingling in the reservoir as a result of replenishment and convection dynamics (Turner and Campbell, 1986; Petrelli et al., 2011; Longo et al., 2012a; Garg et al., 2019).

Arc volcanoes are known for their dominantly explosive character. Their erupted products span the so-called andesitic magmatic suite, ranging from basaltic andesite to rhyolite. The occurrence of repeated pre-eruptive mixing events involving andesitic and dacitic melts is often recognized in the discharged magmas [e.g., (Tamura et al., 2003; Shane et al., 2008; Conway et al., 2020)].

Here we employ the 3D physical model described above to study the physics of mixing between andesite and dacite magmas. We use these magmas because arc-volcanism emits mainly the andesitic suite of magmas. The objective is to study magma dynamics for this suite of magmas through numerical simulations and extract some geologically meaningful information. We remark that the model and numerical scheme employed here is applicable on any suite of magmas. In particular we model the convection dynamics emerging from a gravitationally unstable stratification of andesite and dacite magmas in a shallow reservoir. We refer to a set up that represents a 3D extension of the one employed in previous work (Garg et al., 2019), allowing us to also compare between 2D and 3D dynamics.

We present the results of the numerical simulations by first analysing the 3D convection dynamics, then comparing with the 2D case in (Garg et al., 2019).

As we mention above, the 3D simulation results presented above make use of 2.6 million tethrahedral elements in order to resolve the 3D Rayleigh-Taylor flow structure determined by the adopted initial unstable configuration. Figure 11 shows a zoom view of one isosurface distribution from Figure 2. In particular, superposition of the computational mesh illustrates well the resolution level achieved, and the kind of details that are resolved.

To check the parallel performance of GALES, we conduct a strong scaling test which is widely used to check the ability of a software to deliver results in less time when the amount of resources is increased. The parallel performance is quantified by comparing the actual speedup with the ideal speedup for a given set of processors. The actual speedup and the ideal speedup are defined as Actual speedup = t r t N r < N (10) Ideal speedup = N r (11) where t _r and t _N are the computational times taken by r processors and N processors, respectively. Parallel efficiency is computed as the ratio of the actual speedup and the ideal speedup for a given number of processors: Efficiency = t r t N . r N (12)

An ideal scalable software should result in a linear speedup. However, this is hardly achieved in real situations. For a given mesh the parallel efficiency decreases as we increase the number of processors, mostly as a consequence of increased time spent in communication among processors. The scaling tests are performed on the Marenostrum supercomputer at the Barcelona Supercomputing Centre (BSC) (https://www.bsc.es/marenostrum/marenostrum/technical-information) within the project GALES-3D (2010PA5625) in the frame of a PRACE preparatory access A (https://prace-ri.eu/).

The scaling tests were run on three different meshes with progressively increasing size. The mesh models are listed in Table 2. The simulations use the same setup as described in Section 3.1. The adopted numerical methods and the parallelization strategy employed are illustrated in detail in Appendixes A,B. Here we only recall the monolithic approach to solve the linearized system of equations describing the space-time system dynamics, which is effective in reducing data transfer in parallel computing implementations (Dowd and Charles, 2010).

Figure 12 shows the strong scaling results for each test case. The left panels in the figure display the time taken by the assembly procedure, the linear solver and the total time. We report here the average time taken for a non-linear iteration over 10 time steps (the observed standard deviation being 2 × 10^–2 s). The right panels in the figure plot the speedup and the efficiency as a function of number of cores.

Overall, the scaling behavior of GALES is quite satisfactory in the explored range of computational elements and cores (a parallel efficiency of 0.6 or greater is taken here as satisfactory). The left panels in Figure 12 show that most of the computational time is spent in the assembly of the linearized system of equations. The solution time is only a fraction of the assembly time, increasing in relevance with increasing number of cores thus (right panels) with decreasing parallel efficiency. The total computational time decays nearly linearly (in the log scale of the figures) with increasing number of cores, and the parallel efficiency (right panels) is seen to decrease below satisfactory only for the largest number of cores employed for each different mesh size.

Parallel performance primarily depends on load balancing and inter-processor communication. Our linear solvers are distributed and use peer-to-peer communication to solve the system. Table 3 displays the partition statistics for test case 3 in Table 2, and illustrates well the reasons and conditions under which the parallel efficiency becomes less than satisfactory. Specifically, we display the average number of elements on a processor and the average percentage of inter-processor boundary nodes, computed from the load balanced partitions generated by the Metis software (see the Appendix B for further details). Load balance when increasing the computational resources reduces the overall execution time. However, by increasing the number of processors the percentage of boundary nodes lying at inter-processor boundaries also increases, implying that the processors need to communicate more data. This translates into increased communication time and works towards decreased parallel efficiency. In other words, for any specific problem there is an optimal balance between decreasing load per core and increasing core inter-communication when increasing the overall computational resources. The results in Figure 12 show that for GALES, and in the range of the strong scaling exercise described here, a one order of magnitude decrease in total computational time is allowed before such a balancing condition is achieved.

In this paper we present a 3D parallel code for computing compressible to incompressible multi-component magma dynamics, compare numerical simulations of 3D vs. 2D dynamics for a simple test case represented by the development of Rayleigh-Taylor instability in a stratified, idealized magma chamber, and illustrate code performance by analyzing the results of a strong scaling test involving up to nearly 50 million computational elements and up to > 12,000 computational cores. The computational speedup is close to ideal and above satisfactory levels as long as the element/core ratio is sufficiently large so as to limit boundary nodes to less than half total nodes. The demonstrated parallel efficiency is such as to guarantee efficient use of HPC resources in 3D applications to more realistic magmatic configurations including geometrically complex multiple chamber, dike and conduit systems, likely requiring a number of computational elements exceeding the maximum employed here. Recent extension of GALES to include fluid-structure interaction and coupling with solid elasto-dynamics (Garg et al., 2021) further extends the accessible computational domains requiring exploration of the potential towards exascale computing (e.g., as in the ChEESE European Center of Excellence, www.cheese-coe.eu).

The numerical simulations performed here, designed to compare with previous 2D simulations, illustrate the 3D dynamics of magma chamber overturning due to gravitationally unstable initial conditions and development and evolution of Rayleigh-Taylor instability. The numerical results illustrate to a high resolution level the 3D dynamics of development and growth of the instability, the formation and reciprocal interaction of buoyant plumes of magma, the complexities of the associated vorticity patterns and associated magma mixing, and the evolution towards a stable dynamic state preserving compositional heterogeneities and resulting in a dynamic layering of the magma chamber.

Remarkably, the much simpler 2D simulation approach is found to reproduce the more realistic 3D dynamics on a zero/first-order level, showing similar overall dynamics occurring over comparable time scales, and similar overall effects due to increased magma viscosity. On a more quantitative level, the 3D dynamics show however important differences, and in particular they produce faster plumes leading to more efficient magma mixing than for the corresponding 2D approximation. These results are in general consistent with those from the engineering literature (Young et al., 2001; Cabot, 2006) where similar overall qualitative consistency but important quantitative differences between 2D and 3D simulations are evidenced, although for high Re incompressible flows. In particular, our results suggest that the extent of approximation introduced by the 2D simplification may depend on the specific conditions, and may be smaller for low viscosity conditions. That could be relevant in light of the substantial computational efforts required by the solution of 3D dynamics, and we deserve to evaluate it further, e.g., for more complex geometrical systems over a range of magmatic compositions.