Melt-rock reaction in the asthenospheric mantle: Perspectives from high-order accurate numerical simulations in 2D and 3D

The style and mode of melt migration in the mantle are important to the interpretation of basalts erupted on the surface. Both grain-scale diffuse porous flow and channelized melt migration have been proposed. To better understand the mechanisms and consequences of melt migration in a heterogeneous mantle, we have undertaken a numerical study of reactive dissolution in an upwelling and viscously deformable mantle where solubility of pyroxene increases upwards. Our setup is similar to that described in [1], except we use a larger domain size in 2D and 3D and a new numerical method. To enable efficient simulations in 3D through parallel computing, we developed a high-order accurate numerical method for the magma dynamics problem using discontinuous Galerkin methods and constructed the problem using the numerical library deal.II [2]. Linear stability analyses of the reactive dissolution problem reveal three dynamically distinct regimes [3] and the simulations reported in this study were run in the stable regime and the unstable wave regime where small perturbations in porosity grows periodically. The wave regime is more relevant to melt migration beneath the mid-ocean ridges but computationally more challenging. Extending the 2D simulations in the stable regime in [1] to 3D using various combinations of sustained perturbations in porosity at the base of the upwelling column (which may result from a viened mantle), we show the geometry and distribution of dunite channel and high-porosity melt channels are highly correlated with inflow perturbation through superposition. Strong nonlinear interactions among compaction, dissolution, and upwelling give rise to porosity waves and high-porosity melt channels in the wave regime. These compaction-dissolution waves have well organized but time-dependent structures in the lower part of the simulation domain. High-porosity melt channels nucleate along nodal lines of the porosity waves, growing downwards. The wavelength scales with compaction length and varies slowly and systematically as a function of time. In the presence of a sustained porosity perturbation at the inflow (due to presence of low solidus materials), formation of self-organized wave fields is delayed. A single high-porosity dunite channel first develops, similar to those observed in the stable regime [1]. With increasing time, the amplitude of the wave increases, eventually overwhelms the high-porosity channel. A self-organized wavefield then emerges. The wave pattern depends weakly on the sustained porosity perturbation. The strong correlation in spatial distribution between dunite channel and mantle heterogeneity therefore is lost in the wave regime. In 3D, the idealized dunite channel is cylindrical, not tabular as often inferred from field observations. The latter may be formed in the presence of shear deformation, a subject that we are actively exploring [4]. [1] Schiemenz et al. (2011) Geohpys. J. Int. 186, 641-664. [2]. Bangerth et al. (2007) ACM Trans. Math. Software 33, doi: 10.1145/1268776.1268779. [3] Hesse et al. (2011) Geohpys. J. Int. 187, 1057-1075. [4]. Baltzell et al. (2013) AGU abstract.