Numerical simulation of magma-rock interaction at Krafla volcano using OpenFOAM software and a simplified thermal model

Abstract. We present a 2D numerical modelling study aimed at exploring magma-rock interaction following the emplacement of a magmatic sill into cold shallow crust. An interface-tracking solver was developed, based on the open-source OpenFOAM package that enables simulation of heat and momentum transfer between magmas of different compositions, with contrasting densities, thermal properties, temperatures, crystal contents, and strain-rate dependent viscosities. Two scenarios are considered to reconstruct sharp temperature gradients and explain the presence of fresh rhyolitic fragments excavated from approximately 2 km depth during IDDP-1 drilling at Krafla caldera in 2009: partial melting of felsic crust triggered by either (1) a 300 m thick rhyolite intrusion or (2) a 100 m thick basalt sill. We also assume two possible magma emplacement periods, during the Krafla Fires (1975–1984, ∼ 35 years before drilling) and during the Myvatn Fires (1724–1729, ∼ 300 years before drilling). In Scenario (1), vigorously convective molten rhyolite produces a temperature jump (400 °C) over approximately 25 m (∼ 16 °C m−1) 35 years after emplacement. After 300 years, the thickness of these molten rocks reaches approximately 70 m, however, the thermal gradient becomes too small (less than 5 °C m−1) to explain the IDDP-1 observations. In Scenario (2), because of large density contrasts between the injected basaltic magma and molten rhyolite, two separate convective layers are formed. The thickness of molten rocks reaches about 40 m after 35 years, and the propagating melting front produces a sharp temperature gradient in the undisturbed rocks, greater than 26 °C m−1. These results together with previous petrological studies lead us to conclude that this second scenario of a basaltic intrusion provides a more robust explanation for the extreme geothermal gradient encountered in 2009 than the first scenario. By comparing with a simplified 1D thermal model and performing parametric tests, we argue that both 2D and 1D numerical approaches help constraining better magmatic convection at such extremely high Rayleigh and Prandtl numbers.