High-Performance Computing of 3D Magma Dynamics, and Comparison With 2D Simulation Results

The dynamics of magma is often studied through 2D numerical simulations because 3D simulations are usually complex and computationally expensive. However, magmatic systems and physical processes are 3D and approximating them in 2D requires an evaluation of the information which is lost under different conditions. This work presents a physical and numerical model for 3D magma convection dynamics. The model is applied to study the dynamics of magma convection and mixing between andesitic and dacitic magmas. The 3D simulation results are compared with corresponding 2D simulations. We also provide details on the numerical scheme and its parallel implementation in C++ for high-performance computing. The performance of the numerical code is evaluated through strong scaling exercises involving up to > 12,000 cores.


INTRODUCTION
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 liquiddominated 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 twodimensional, 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.

MAGMA FLOW MODEL
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 multicomponent 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 fluidstructure 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 2 O and CO 2 . The mixture is assumed to be in chemical, mechanical and thermodynamic equilibrium. The flow is governed by the following equations: 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 threedimensional 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: 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 /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 where D is the mass diffusion coefficient matrix. Viscous flux is modelled as 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: 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 0 v 0 /l 0 , where t 0 is the relaxation time of the crystal, v 0 is the flow velocity and l 0 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 14 m −3 (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 spacetime 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 2 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 2 O + CO 2 ) gasmelt equilibrium modelling (Papale et al., 2006). The viscosity of each magma component is modelled as a function of oxide composition, dissolved H 2 O and temperature (non-Arrhenian) as described in (Giordano et al., 2008), with the effect of nondeformable 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 − α 2 T/(βρ), where α and β are isobaric expansion coefficient and isothermal compressibility coefficient, respectively. In this work we use a constant thermal conductivity, κ 1.2 Wm −1 K −1 (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.

APPLICATIONS TO MAGMA MIXING DYNAMICS
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.

Simulation Setup
The computational domain represents a prolate ellipsoid with longer axis in the x direction ( Figure 1). The centre of the ellipsoid is placed at 4.1 km depth. The lengths of the semiaxes in x, y, and z directions are 400, 100 and 100 m respectively. With respect to the 2D setup in (Garg et al., 2019), the third dimension added here (z-axis) is short enough to cause significant deviation from 2D conditions (approximated when the neglected dimension is much longer than the considered ones). We perform two simulations which correspond to run cases #1 and #3 in (Garg et al., 2019). For consistency with the corresponding 2D cases in (Garg et al., 2019), we keep the same numbering throughout in this work ( Table 1). In the simulations, andesite at a temperature of 927°C is placed at the bottom of the domain, while the upper part is filled with dacite at a temperature of 876°C. A horizontal interface, separating the two magmas, is set at 4,150 m depth (Figure 1). The chemical composition of the two magmas and their volatile contents are reported in Table 1. The two cases differ for only magma viscosity, with case #1 corresponding to locally defined, spacetime dependent viscosity computed as described above, and case #3 equal to case #1 with the viscosity arbitrarily multiplied by a factor 10 everywhere in space and time. The effects of varying the amount of volatile contents in the 2D setup were studied in (Garg et al., 2019). While such a viscosity increase may approximate the effect of non-reactive crystals, which affect viscosity much more than other properties including density, the aim here is mostly that of evaluating the 3D dynamics and comparing with the 2D case, over a range of viscosities thus of dynamic time scales; as well as that of evaluating the 3D code performance in terms of computational time and scalability. The initial pressure distribution is computed by considering the lithostatic load at the chamber roof (average rock density 2,500 kg/m 3 ) and a horizontally uniform magmastatic profile. The bottom panel of Figure 1 (red lines) displays the initial pressure and density profiles along the chamber axis. No-slip adiabatic conditions are employed at the chamber walls. The numerical scheme presented in the Appendix A is employed with 2.6 million linear tetrahedral elements. The side length of the computational elements ranges 1-4 m. The numerical integrals are computed with four Gauss quadrature points. The unstructured tetrahedral mesh results in uneven interface providing the initial perturbations that destabilize the interface between the two magma types. The simulations are run on an HP cluster system at INGV Pisa, composed of 432 Intel Xeon 2.3 GHz cores distributed over six nodes connected through a low latency 100 Gbps Infiniband network. Scaling tests reported below and involving up to thousands cores for a much shorter computational time are run instead on the supercomputing facilities at the Barcelona Supercomputing Center.

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

3D Dynamics
The 3D dynamics are illustrated here mostly through comparison between the two simulation cases 1 and 3 in Table 1, with the latter being identical to the former, except for the computed viscosity which is everywhere in space and time arbitrarily multiplied by a factor 10. We anticipate that as for the 2D case (Garg et al., 2019), the more viscous situation corresponds to slower convection dynamics, lower number of buoyant plumes of andesite-rich magma rising through the dacite, and finally lower mixing efficiency. Figure 2 illustrates well such differences. The figure shows the evolution of the isosurface corresponding to a mass fraction of andesite equal to 0.5 (which at time zero corresponds to the interface between the two magma types) for cases 1 and 3. The differences are striking: after 100 s case 1 shows a highly dynamic state with a complex structure made of several inter-digitized rising plumes interacting with each other and occupying the entire chamber; while at the same time, only minor perturbations appear on the interface initially separating the two dacitic and andesitic magmas. At a later time when the dynamics are well developed also for case 3, this case displays a much simpler overall structure with many less plumes, each one on average much bigger than for case 1. The 0.5 mass fraction isosurface for the more viscous case 3 is conserved within the chamber over a much longer time, and it is still visible after more than 2 h from start of the simulation. On the contrary, the same isosurface is completely lost in case 1 after only 5 min, as a consequence of much faster and more intense magma mixing for this lower viscosity case. Figure 3 highlights the differences in plume structure for the two cases, through a planar view of the interface at an early stage of its destabilization. The higher the viscosity (case 3), Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 760773 the lower the number density (and the larger the average size) of the rising plumes formed as a consequence of Rayleigh-Taylor Instabilities. Figure 4 shows instead the distribution of velocity in the uprising plumes, at the two times (90 and 250 s, respectively) at which the maximum velocity is attained for the two cases 1 and 3. That maximum velocity corresponds to 8.2 m/s in case 1, while it is only 3.5 m/s for the more viscous case 3.   Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 760773 6 Figure 5 illustrates the evolution of mixing throughout the computational domain. Case 1 corresponds to viscosity computed by model described in (Giordano et al., 2008) and case 3 is with tenfold increase in viscosity than case 1. Initially, for both cases 1 and 3 the computational nodes host either pure dacite or pure andesite. In both cases the less viscous, less abundant andesitic end-member quickly vanishes as a pure component (at the scale of the resolution of the present simulations, which is of order 1 m), whereas nearly pure dacite continues to be largely present in the system at the latest simulation times. Faster and more efficient mixing for case 1 is clearly visible as an earlier decrease of the maximum length of the compositional bars (that is, faster decrease of the number of nodes hosting pure dacite) as well as narrower compositional interval span inside the chamber at latest simulation times, when efficient convection is terminated in both cases (further illustrated below). For the more viscous case 3 the andesite can be seen to constitute at most about 50% of individual computational nodes, whereas for case 1 the maximum proportion of andesite in individual computational nodes at process end is only about 30%.
In both cases 1 and 3 by the latest times, the system is decompressed throughout by 6 and 2 bars, respectively, and the density evolves from an initial step profile to a smooth one (Figure 1). For both simulation cases 1 and 3, pressure and temperature distributions along z 0 and x 0 planes is provided in the Supplementary Material. We also display the profiles of scaled temperature and composition along the chamber axis in the Supplementary Material.

Comparison Between 3D and 2D Simulation Results
As it is explained above, the present 3D simulation cases correspond to previous 2D simulations in (Garg et al., 2019), so to confidently explore the effects of 3D vs. 2D simulation setup. In particular, the numerical code and the physical and numerical setup in the 2D and 3D cases are the same, except for the specific aspects defining 2D vs. 3D simulated dynamics. The average length of computational elements across the initial compositional interface is of order 1 m in both cases, and the number of elements in the 2D cases and along domain-centered 2D slides in the 3D cases are of order 10 5 in both cases.
One of the major results from the simulation of 2D magma convection dynamics was the achievement of a stable dynamic state characterized by separate convective regions effectively impeding further mechanical mixing (Garg et al., 2019). Figure 6 shows that such a major result is confirmed by 3D numerical simulations: dynamically stable, largely closed circulation patterns separating regions with different composition emerge, explaining the long-term maintenance of the heterogeneities at advanced times in Figure 5. In fact, the circulation patterns in Figure 6 give rise to a stable, dynamic layering inside the magma chamber (Figure 7). The details of the circulation patterns are not easy to compare, as those patterns are intimately three-dimensional for the 3D simulations, therefore, comparing a slice cut with the 2D case, e.g. as in Figures 8, 9 below, would not make any sense. Here we stress the overall first order agreement between the 2D and 3D simulation cases with respect to the major conclusion that magma chamber overturning and associated convection and mixing dynamics lead to the achievement of a stable dynamic state preserving compositional heterogeneities over the long term. We discuss in (Garg et al., 2019) some major consequences for the interpretation of compositional heterogeneities in magmas erupted during individual eruptions. A more direct comparison between the simulated 2D and 3D dynamics is displayed in Figures 8, 9 for the simulation cases 1 and 3, respectively. For the 2D case each panel in the figures represents the entire computational domain, with the third neglected dimension (perpendicular to the sheet) assumed much longer than the two simulated ones so as to satisfy the 2D approximation. On the contrary, for the 3D case each panel reports a vertical slice cut across the center of the 3D computational domain, with the third dimension, simulated but not visible from the slice cut, being equal in size to the vertical dimension (Figure 1).
At first sight, the 3D and 2D dynamics in Figures 8, 9 appear quite similar, especially in light of the different evolutions characterizing the less and more viscous cases, respectively, one and 3. That may appear surprising, considering the small length of the third chamber dimension in the present 3D simulations. However, such qualitative similarities are consistent with previous findings, e.g., (Young et al., 2001;Cabot, 2006) (although those authors investigate high-Re  incompressible flows), and as for those cases, we show here that the differences are also relevant. The comparison shows that the interface destabilization time in the two 3D and 2D cases is very similar [while it is significantly influenced, as any other aspect of the overall dynamics, by magma viscosity. The role of viscosity-and of volatile contents-is however described in (Garg et al., 2019), and it is only marginally considered here]. However, the 3D geometry of the interface appears to result in more complex perturbation structure comprising a broad range of scales, compared to the geometrically much simpler 2D perturbations. At later times such a richer 3D structure is visible as a much less symmetric geometry  Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 760773 9 of the rising plumes, and more widespread plume size range. The plumes in the 3D cases also appear to rise faster, and are more mixed than for the 2D cases. That is well evident at time 100 s for case 1 (Figure 8) or time 200 s for case 2 (Figure 9): the 2D plumes maintain a much larger proportion of nearly pure andesite (dark red), which is instead much less represented (Figure 9), or practically absent (Figure 8), in the 3D cases at the same time. At the latest simulation times corresponding to achievement of a stable dynamic stratification as discussed above, the more viscous case in Figure 9 shows that andesite-rich (70-80 wt%) magma concentrates close to magma chamber top in the 2D case, whereas only < 60 wt% andesitic magma occupies the same chamber region in the 3D case.
The faster, more efficient mixing dynamics found for the 3D case are best highlighted in Figure 10. Here, to establish a measure of mixing, we refer to progressive reduction of the overall compositional heterogeneity inside the chamber. As a quantitative measure of such heterogeneity we compute the standard deviation (σ) of the mass fraction of andesite in the overall computational domain. For any time we first determine the mean value of the composition in the entire domain, then use that value to compute (σ), which measures the extent of dispersion of the composition around its mean value. The larger the value of (σ), the lower the extent of mixing: where N is the total number of nodes in the computational domain, x is mass fraction, and the horizontal bar indicates the mean value over the entire computational domain.
The computed evolutions of σ in Figure 10 for the 2D and 3D simulations evidence the different stages of the overturning process (Garg et al., 2019), starting in all cases with a low slope section corresponding to the initial phase 1 of development of the Rayleigh-Taylor instability, followed by a high slope section corresponding to highly efficient mixing during the convective phase of plume rise and vortex formation (phase 2), and terminating with an essentially flat section (phase 3) when the final stable dynamic state described above, preventing significant further mixing, is achieved. Faster development of the Rayleigh-Taylor instability and transition to phase 2 (efficient plume rise) is highlighted in the zoom view of Figure 10B. Similarly, faster overall dynamics and more efficient mixing for the 3D case translate in earlier achievement of the flat section corresponding to phase 3. Note that for the high-viscosity case 3, the 2D simulation never attains a completely flat section, indicating that some mixing continues to be effective up to the last simulated time > 2 h. That suggests that constraining the magma circulation patterns over a 2D plane (or better, assuming zero gradient of any flow variable including velocity along the neglected dimension as it is implicit in 2D simulations) may cause the flow streamlines to distort to an extent sufficient to break the closed circulation patterns in Figure 6, resulting in further mixing not seen at such late times from the 3D simulations. Finally, and mostly relevant, the final extent of homogenization (that is, the overall change in σ) is significantly larger for the 3D simulation cases, reflecting enhanced convection and mixing efficiency with respect to the 2D case. The difference is important, amounting to 18% of the total σ change computed from the 2D simulation for case 1, and as large as 55% of that change for case 3. Therefore, not only mixing and homogenization in the magma chamber are significantly enhanced when considering the flow dynamics in a more realistic 3D environment; also, the extent of such enhancement depends substantially on the specific conditions. While the less viscous case 1 leads to more homogenization, the extent of change when moving from 2D to 3D turns out to be larger for the more viscous case 3. In other words, it seems plausible, based on our first 3D simulations and comparison with the corresponding 2D cases, to suggest that the errors and approximations introduced by neglecting more realistic 3D dynamics may vary substantially depending on the conditions considered, and may increase in relevance with increasing magma viscosity.

CODE PERFORMANCE
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 FIGURE 11 | Zoom view of the computational mesh and isosurfaces Y 0.5 at t 100 s for case #1.
Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 760773 given set of processors. The actual speedup and the ideal speedup are defined as Actual speedup t r t N r < N (10) 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: An ideal scalable software should result in a linear speedup. However, this is hardly achieved in real situations. For a given  1  175001  1004750  1050006  96  96-1,536  2  1951658  11782329  11709948  192  192-6,144  3  7802061  47581941  46812366  384  384-12,288 FIGURE 12 | Strong scaling results for test cases in Table 2.
Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 760773 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.

DISCUSSION AND CONCLUSION
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 elastodynamics (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.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
DG and PP designed the study, analysed and discussed the results, and wrote the article, DG carried out the numerical simulations and postprocessed the results.

ACKNOWLEDGMENTS
We gratefully acknowledge the reviewers for their comments, which improve the quality of the manuscript. We acknowledge that the results on the code performance have been achieved using the PRACE Research Infrastructure resource MareNostrum based in Spain at Barcelona Supercomputing Centre (BSC). One of us (D.G.) was supported by an EPOS-IT grant. Conflict of Interest: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

APPENDIX A NUMERICAL METHOD
The system of Eqs (1-4) can be written and solved in a fully coupled monolithic manner as a single transport system (Hauke and Hughes, 1998;Longo et al., 2012b;Garg et al., 2018a): where the vectors are given by In the above equations the following notions are used: δ i δe i and τ i τe i , where e i is the unit vector in the i th direction and δ [δ ij ] is the Kronecker delta. The sub-indexes i and j stand for spatial coordinates, i.e. i, j 1,2,3; while k stands for the mixture component. For spatial coordinates the Einstein summation convention of repeated indexes is used throughout.
Equation (13) can be rewritten for any independent set of variables X (Shakib et al., 1991;Hauke and Hughes, 1994) as where A 0 U ,X , A i F a i,X is the ith Euler Jacobian matrix and K [K ij ] is the diffusivity matrix with K ij X ,j F d i . The monolithic approach used in this work has several advantages over segregated approaches in terms of simplicity of the formulation, robustness of the solution approach, and smaller data transfer when doing parallel programming. In the present study, the vector X has been chosen as the set of pressure primitive variables, i.e. X [p, v, T, Y k ], which allows modelling of both compressible and incompressible flows by a unified formulation (Hauke and Hughes, 1994;Hauke and Hughes, 1998;Longo et al., 2012b;Xu et al., 2017;Garg et al., 2018a;Garg et al., 2018b;Garg et al., 2021).
The boundary value problem can be expressed as below. Consider an open spatial domain Ω with boundary Γ, such that Γ Γ G ∪ Γ H and Γ G ∩ Γ H ϕ, where Γ G and Γ H are the Dirichlet and Neumann parts of the boundary, respectively. The strong form of the problem consists of finding the solution vector X : Ω → R neq , where n eq is the number of equations of the system, such that for the given essential boundary conditions X G and the natural boundary conditions X H , the following equations are satisfied: where R(X) is the residual of the equations and L represents the transient-advective-diffusive operator such that The weak form of the above equations can be expressed as: given a trial function space T {X | X ∈ (H 1 ) neq , X X G on Γ G } and weighting function space which by substituting the definition of R(X), integrating by parts and applying the boundary conditions, can be written as

A.1 Variational Multiscale Formulation
We consider the transport operator L in Eq. 17 as quasi-linear and solve Eq. 13 by variational multiscale method (VMS) (Hughes, 1995;Hughes et al., 1998;Bazilevs et al., 2007). The VMS method has been successfully used for compressible flows (Franco and Rafael Saavedra, 2006;Rispoli et al., 2015) and turbulent flows (Bazilevs et al., 2007) without employing any ad hoc terms such as eddy viscosities. In the VMS method the solution vector X is decomposed into the resolved (coarse, gridscale) finite element solution X and unresolved (fine, sub-grid) error X9, Similarly W is decomposed as The VMS method consists of substituting the above splitting into the weak form and solving the following two subproblems for coarse-scale and fine-scale: The essence of VMS method is to solve for the coarse scale solution numerically and compute the fine scale part either analytically or approximate it through an algebraic expression. With this aim, using adjoint duality, Eq. 22 can be written as Assuming X9 0 on element boundaries, Eq. 23 can be expressed as: The solution of the above equation can be represented as a function f of X and R( X): Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 760773 X′ f X, R X To obtain a simple, basic and computationally efficient method, the VMS formulation relies on approximating X9 with the product of element-wise algebraic stabilisation operator P and the coarse scale residual, R( X) (Bazilevs et al., 2007), Complete VMS formulation along with the discontinuity capturing operator (see below) (Tezduyar and Senga, 2006) can be expressed as:

A.2 Stabilization Operator
In the present work we follow the same design for the P matrix as it was developed in (Hauke, 2001) and extended to multicomponent conditions in (Longo et al., 2012b). We first design P U for the conservation variables and then transform it into the pressure-primitive-variable formulation (P X ) using the following expression P P X X ,U P U (30) P U is given as: P U diag P c , P m , P m , P m , P E , P Y k where the diagonal entries are given by where Δt is the time step, h e is the element size along the streamline direction, c is the local sound speed and d k is the mass diffusion coefficient of kth component.
To handle incompressibility the entry of the first row of the P matrix is modified as P c P −1 c + ρP m * g · g −1 −1 where P c and P m are the diagonal entries of P matrix corresponding to the continuity and the momentum equations, and