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

Gerbault, Muriel; Melnik, Oleg; Borisova, Anastassia

doi:https://doi.org/10.5194/egusphere-2025-3036

Preprints

https://doi.org/10.5194/egusphere-2025-3036

Preprints

07 Jul 2025

| 07 Jul 2025

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

Muriel Gerbault, Oleg Melnik, and Anastassia Borisova

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 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 meters (~16 °C/m) 35 years after emplacement. After 300 years, the thickness of these molten rocks reaches approximately 75 m, however, the thermal gradient becomes too small (less than 5 °C/m) 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 exceeds 30 m after 30 years. The rapid melting front propagation causes a sharp temperature gradient in the undisturbed rocks (28 °C/m). We conclude that the second scenario provides a more reliable explanation for the existing data and is well supported by previous petrological studies. By comparing with a simplified 1D thermal model and performing parametric tests, we argue that our numerical approach is suitable for studying magmatic convection at such extremely high Rayleigh and Prandtl numbers.

Received: 25 Jun 2025 – Discussion started: 07 Jul 2025

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.

Download & links

Muriel Gerbault, Oleg Melnik, and Anastassia Borisova

Status: final response (author comments only)

RC1:
'Comment on egusphere-2025-3036', Catherine Annen, 25 Jul 2025

General comment
The study presented in this manuscript seeks to account for the high temperature gradient revealed by IDDP-1 drilling in Krafla caldera. The study consists in a series of numerical models involving a basalt or rhyolite sill that is convecting and melts the felsic crust above. A best fit was found with observation with a 100m sill basalt although numerical issues were encountered depending on the chosen viscosities.
I like that the authors report both successful and failed models as negative results are also important to the advance of science. The results are interesting and coherent. The text is not always clear, but it is often more an issue with the style than with the science. I found the discussion particularly difficult to follow: it lacks a clear narrative and some sentences are non sequitur.
Specific comments
I have no major issues with the science as the model and its limitations are well explained.
The only point is that although I understand the use of the Myvatn Fires and Krafla fires to have some timescales for the model, I am not sure that there is any strong argument for the timing of the sill emplacement to coincide with an eruption. It might be the opposite with a sill emplacement being a failed attempt of the magma to reach the surface.
Detailed comments and corrections
l.35-37: Couldn’t the high temperature gradient also be explained by a very recent intrusion and heat that hadn’t had time to diffuse?
l.44-45: I suggest removing “Until today”. Can you be more specific as why a single reservoir cannot explain bimodality (just one or two sentence)? Does it mean that the rhyolite has to be melted crust but if bimodal reservoir was possible then the rhyolite could have come from the reservoir?
l.49: can the crust be felsite? I thought felsite was for volcanic rock. Maybe “rock” would be better than crust here.
l.54 It is not clear what is the process advocated by Simakin and Bindeman. If it is precisely what has been described before, I suggest giving the reference and removing “corresponding to the process advocated by”
l.95-100: What is the dimension of C_i?
In table 1, crystal densities are lower than melt densities, isn’t it the other way round?
l.130: I find confusing to have k for thermal diffusivity and κ for thermal conductivity as it is usually the other way round.
l.163: there is a contradiction between text and figure regarding the bottom boundary condition: no slip vs free slip.
Equation (5): The meaning of C_i is not clear to me and why the sum is 1 to 3, when above on l.96, i is 1 or 2. Is C_i in equation (8) the same? It seems to be a temperature there.
l.197: add “a”: injected into a cold rhyolite crust.
l.206: “increasing” instead of increases.
l.206: rephrase: a jump cannot form a plateau. Maybe something like: “A plateau in temperatures follows an initial sharp thermal jump at the basalt-rhyolite boundary. The plateau corresponds to the convective layer of partially molten rhyolite.”
l.209 and 211: the term “evolves” is more suitable for a time variation than a spatial variation. A possible alternative phrasing could be: “After 5 years, the temperature across molten/unmolten rhyolites jumps from about 500 to 1000oC over a depth [..]”
l.215 I suggest removing “Naturally”
l.217 Instead of “tendency can be draw”, “tendency can be identified”.
l. 217: Needs rephrasing. If I understood correctly (not sure it is the case): “the difference in the melting zone thickness reaches 25% for a cell size change from 0.5 m to 1 m (MZT 20 m vs 27 m); it reduces to 15% for cell sizes change from 0.25 to 0.5 m, etc.
l. 223: I suggest replacing “Therefore one can expect it would actually reach [..] if we would reach the critical thickness” with “Therefore, we hypothesize that it would reach [..] if the critical mesh size could be reached”
l. 226: “that” instead of “how”
l.260: the sentence is ambiguous. Is it 2 TBLs, one in each magma, or 4 TBLs, two in each magma?
l.263: “same”: what is same referring to?
l.283: “this would also affect the melting front thickness by a comparable amount.” I don’t understand this sentence; comparable to what?
l.284-285: I suggest rephrasing: “The intrusion’s thickness controls the amount of melting of the overlying rhyolite but is not well constrained by geophysical data.”
l.290: The sentence sounds awkward, consider rephrasing.
l.294: What does “greater modelling time mean” (longer, shorter)?
l.295-296: Sentence starting with “Beside” is not clear.
l.299-300: Awkward phrasing. A choice cannot be sufficient.
l.302: “the water content of Krafla’s fresh lavas to granophyres is low” I don’t undersand. Do you mean fresh lavas and granophyres? or are fresh lavas granophyres?
l.303: what does subscript VSMOW mean in relation with ¹⁸O. Please, clarify for non-experts.
l.303-304. Sentence stating with “Thus” I do not understand what is meant. It sounds non-sequitur.
l.333-335: This sentence is difficult to follow and lack specificity: segregation of pure partial melt from where to where? What is meant by pure? Assimilation of what? Fractional crystallization of what?
l.336: What is the Víti felsite?
l.343: It is not clear how the last sentence of the section (starting with “In this context”) links to the former statements (d¹⁸O heterogeneities).
l.351-352: I am not convinced that the sill intrusion has to coincide with an eruption.
l.465: What do you mean by one can set a scaling factor. How does it link with eq. A2. Please, be more specific.
l.482: What is the curve shape? What curve?
l.498: “it can destabilize gravitationally” Is it really unrealistic? Do you mean unrealistic in this specific case or in general? It looks like stoping.

Caption of figure 1: the geophysical anomaly is said to be purple but I see it blue.
Caption of figure 2: Have T_inf the same in caption and figure. Delete “hot”? a) Case of rhyolite -ADD/ with composition- similar to [..]
Figure 4: I am confused by the diagrams of velocity. The velocity seems to be non zero in areas that are solid. How is that?
Figure 5 and 8 top: why different colours for the symbols?
Figure 5 middle: We lose the value of mean velocity after 400 yrs. Would you consider changing the scale of Y axis?
Caption figure B3: The last sentence sounds too colloquial.

Citation: https://doi.org/10.5194/egusphere-2025-3036-RC1
- CC1:
  'AC reply on RC1', Oleg Melnik, 13 Aug 2025
  
  Dear Catherine,
  
  Thank you for careful reading and high ranking of our manuscript. Below, we reply to some of your specific comments. We agree with most of the minor comments and will correct the manuscript accordingly.
  The only point is that although I understand the use of the Myvatn Fires and Krafla fires to have some timescales for the model, I am not sure that there is any strong argument for the timing of the sill emplacement to coincide with an eruption. It might be the opposite with a sill emplacement being a failed attempt of the magma to reach the surface.
  We think that an injection of a basaltic or rhyolitic sill after Krafla Fires is highly unlikely because the volcano was very well monitored, and such event will be definitely recorded. We've failed to find any evidence of the post Krafla Fires emplacement. The case of Myvatn Fires is less constraint. There could have been undocumented failed eruptions. We will add comments on this subject to the manuscript.
  l.35-37: Couldn’t the high temperature gradient also be explained by a very recent intrusion and heat that hadn’t had time to diffuse?
  
  That is related to the previous comment. We believe that at list after Krafla Fires the emplacement is not likely to occur.
  ll.44-45: I suggest removing “Until today”. Can you be more specific as why a single reservoir cannot explain bimodality (just one or two sentence)? Does it mean that the rhyolite has to be melted crust but if bimodal reservoir was possible then the rhyolite could have come from the reservoir?
  A recent injection of the rhyolite will not produce enough heat to melt the crust above the sill unless the thickness of a rhyolite sill is significant. Such sill should be detected by geophysical methods. It is not possible to sustain high temperature gradient without active melting of the crust.
  l.95-100: What is the dimension of Ci?
  
  Ci is the concentration of magmas in the VOF method, dimensionless.
  In table 1, crystal densities are lower than melt densities, isn’t it the other way round?
  
  Will be corrected.
  l.130: I find confusing to have k for thermal diffusivity and κ for thermal conductivity as it is usually the other way round.
  
  Will be corrected.
  
  l.163: there is a contradiction between text and figure regarding the bottom boundary condition: no slip vs free slip.
  
  No slip
  Equation (5): The meaning of Ci is not clear to me and why the sum is 1 to 3, when above on l.96, i is 1 or 2. Is Ci in equation (8) the same? It seems to be a temperature there.
  
  Should be “2”
  l.465: What do you mean by one can set a scaling factor. How does it link with eq. A2. Please, be more specific.
  
  We have done two viscosity corrections to see the influence of this parameter. First a multiplication of the viscosity by the scaling factor, second - setting the upper and the lower bounds of the viscosity.
  l.498: “it can destabilize gravitationally” Is it really unrealistic? Do you mean unrealistic in this specific case or in general? It looks like stoping.
  
  It depends on the timescale. If we model only 30-300 years of the system evolution, than the viscosity of host rocks can be set to an unrealistically low value (~ 10^10 Pa s) because no significant viscous deformation occurs in the cold region of the domain. If the viscosity is ~ 10^9 Pa s, diapires develops in the cold rocks which is not realistic on these timescales. Longer timescales will require higher viscosity threshold.
  Figure 4: I am confused by the diagrams of velocity. The velocity seems to be non zero in areas that are solid. How is that?
  
  Velocities of 10^-10 m/s results in total displacement of ~1 m in 300 years. This is negligible in comparison with the domain size and displacements in the convective part of the domain.
  Best regards,
  
  Muriel, Oleg and Anastassia
  
  Citation: https://doi.org/10.5194/egusphere-2025-3036-CC1
  - AC2: 'Reply on CC1', Oleg Melnik, 10 Sep 2025
    
    Dear Catherine,
    
    Thank you for careful reading and high ranking of our manuscript. Below, we reply to some of your specific comments. We agree with most of the minor comments and will correct the manuscript accordingly.
    
    The only point is that although I understand the use of the Myvatn Fires and Krafla fires to have some timescales for the model, I am not sure that there is any strong argument for the timing of the sill emplacement to coincide with an eruption. It might be the opposite with a sill emplacement being a failed attempt of the magma to reach the surface.
    
    We think that an injection of a basaltic or rhyolitic sill after Krafla Fires is highly unlikely because the volcano was very well monitored, and such event will be definitely recorded. We've failed to find any evidence of the post Krafla Fires emplacement. The case of Myvatn Fires is less constraint. There could have been undocumented failed eruptions. We will add comments on this subject to the manuscript.
    
    l.35-37: Couldn’t the high temperature gradient also be explained by a very recent intrusion and heat that hadn’t had time to diffuse?
    
    That is related to the previous comment. We believe that at list after Krafla Fires the emplacement is not likely to occur.
    
    ll.44-45: I suggest removing “Until today”. Can you be more specific as why a single reservoir cannot explain bimodality (just one or two sentence)? Does it mean that the rhyolite has to be melted crust but if bimodal reservoir was possible then the rhyolite could have come from the reservoir?
    
    A recent injection of the rhyolite will not produce enough heat to melt the crust above the sill unless the thickness of a rhyolite sill is significant. Such sill should be detected by geophysical methods. It is not possible to sustain high temperature gradient without active melting of the crust.
    
    l.95-100: What is the dimension of Ci?
    
    Ci is the concentration of magmas in the VOF method, dimensionless.
    
    In table 1, crystal densities are lower than melt densities, isn’t it the other way round?
    
    Will be corrected.
    
    l.130: I find confusing to have k for thermal diffusivity and κ for thermal conductivity as it is usually the other way round.
    
    Will be corrected.
    
    l.163: there is a contradiction between text and figure regarding the bottom boundary condition: no slip vs free slip.
    
    No slip
    
    Equation (5): The meaning of Ci is not clear to me and why the sum is 1 to 3, when above on l.96, i is 1 or 2. Is Ci in equation (8) the same? It seems to be a temperature there.
    
    Should be “2”
    
    l.465: What do you mean by one can set a scaling factor. How does it link with eq. A2. Please, be more specific.
    
    We have done two viscosity corrections to see the influence of this parameter. First a multiplication of the viscosity by the scaling factor, second - setting the upper and the lower bounds of the viscosity.
    
    l.498: “it can destabilize gravitationally” Is it really unrealistic? Do you mean unrealistic in this specific case or in general? It looks like stoping.
    
    It depends on the timescale. If we model only 30-300 years of the system evolution, than the viscosity of host rocks can be set to an unrealistically low value (~ 10^10 Pa s) because no significant viscous deformation occurs in the cold region of the domain. If the viscosity is ~ 10^9 Pa s, diapires develops in the cold rocks which is not realistic on these timescales. Longer timescales will require higher viscosity threshold.
    
    Figure 4: I am confused by the diagrams of velocity. The velocity seems to be non zero in areas that are solid. How is that?
    
    Velocities of 10^-10 m/s results in total displacement of ~1 m in 300 years. This is negligible in comparison with the domain size and displacements in the convective part of the domain.
    
    Best regards,
    
    Muriel, Oleg and Anastassia
    
    Reply
    
    Citation: https://doi.org/10.5194/egusphere-2025-3036-CC1
    
    Citation: https://doi.org/10.5194/egusphere-2025-3036-AC2
RC2:
'Comment on egusphere-2025-3036', Alain Burgisser, 20 Aug 2025

The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2025/egusphere-2025-3036/egusphere-2025-3036-RC2-supplement.pdf

Citation: https://doi.org/10.5194/egusphere-2025-3036-RC2
- AC1: 'Reply on RC2', Oleg Melnik, 01 Sep 2025
  
  Thank you very much, Alain, for this positive and motivating review. We answer below the modifications we made in order to clarify the points you raise, and we hope now that our work is more convincing.
  We will address all minor comments in the updated manuscript, and here only reply to a few points that require clarification from our side.
  186 I do not understand why the grid size selection is evaluated against the melt zone thickness, whereas the Section 2.3 states on l. 155 that the best measure is the Nusselt number. Please clarify why the finding of Stevens et al. (2013) is set aside here.
  We use melting zone thickness for the comparison of different melt resolutions as it is a key outcome of the paper - how much molten rhyolite is produced. Nusselt number gives only first order estimates for required melt size because the parametrization of Stevens et al. (2013) was done for slightly different setup - constant viscosity, one fluid, thermal convection only. Thus we consider that independence of meting efficiency on the melt size is a more important for this study.
  214 (and also l. 290) It is important to state the reason(s) behind the crash: vanishing time step, stalling of the residuals, temperature runaway, instability of Eq (4), ... Currently, I can only infer that the crash probably does not occur because of too high a viscosity (see next comment).
  Indeed, the time-steps diminish during the first ca. 10 years of the modelled run-time, and in successful cases it manages to rise back up again once the convective layers are stably growing. In some cases, the time-step decreases to a very small value and never rises back up again. In other cases in turn, the time-step would re-increase so sharply that it is at some point forced to reduce back down sharply and again, it would get stuck into not being able to re-increase. Sometimes bounding the time-step upper limit helps to overcome the problem. Higher mesh sizes tended to produce more stable runs despite their longer total duration, indicating that the resolution of the physical instability is achieved. At later stages when basalt becomes viscous and the viscosity contrast between the magmas becomes small the mixing intensity increases and, of course, this process is not captured by the VOF method.
  “although the 1D model shows temperature gradients within convective regions”. This is a major drawback as the variable of interest is the temperature gradient. Just looking at Fig. 9, I guess that the 1D gradients across the melt zone thickness are basically meaningless because they are so far away from the 2D gradients. If my guess is correct, I suggest adding a few sentences about this issue as it shows clearly the added value of this 2D study.
  We agree that simplified 1D model does not capture all complexity of the 2D simulations. The important message is that here we directly parameterize the effective thermal conductivity and obtain high values of Nu number. The presence of significant temperature gradients in melt zones mean that 1D effective conductivity (order of 200 W/m/K) is smaller than the effective conductivity in 2D and 2D model resolves the heat transfer on required time and length scales. Temperature gradients in unmolten part above the rhyolite layer and its thickness are well captured by the 1D model.
  330-344. The petrology § are not linked to the framework of this study. Please ensure that it is the case. For instance, fractional crystallization is incompatible with the model assumptions, which needs to be mentioned. Another example is that hydrothermal fluids were needed to produce the partial rhyolitic melts, but 1) the model ignores hydrothermal fluids and 2) the whole § on l. 301-309 is dedicated to show that hydrothermal fluids played no direct role in [...] partial melting and the following reaction of the rock with [putative] basaltic magma.
  We will rewrite this section to cover aspects that are compatible with the current model and show model limitations in comparison with the petrological complexity.
  Best
  Muriel, Oleg and Anastassia
  
  Citation: https://doi.org/10.5194/egusphere-2025-3036-AC1
EC1: 'Comment on egusphere-2025-3036', Virginie Pinel, 08 Sep 2025

Dear authors, Two reviewers have recognised the great interest and added value of your article submitted to Solid Earth and have provided comments and suggestions for modifications. You have already begun to respond to these in the discussion. Your article requires some minor revisions before it can be accepted. I encourage you to submit a revised version (with the changes highlighted) and a detailed response to the two reviews by C. Annen and A. Burgisser.

Best regards

Citation: https://doi.org/10.5194/egusphere-2025-3036-EC1

Muriel Gerbault, Oleg Melnik, and Anastassia Borisova

Viewed

Total article views: 725 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	BibTeX	EndNote
635	74	16	725	44	44

HTML: 635
PDF: 74
XML: 16
Total: 725
BibTeX: 44
EndNote: 44

Views and downloads (calculated since 07 Jul 2025)

Month	HTML	PDF	XML	Total
Jul 2025	176	28	7	211
Aug 2025	269	32	4	305
Sep 2025	190	14	5	209

Cumulative views and downloads (calculated since 07 Jul 2025)

Month	HTML	PDF	XML	Total
Jul 2025	176	28	7	211
Aug 2025	269	32	4	305
Sep 2025	190	14	5	209

Viewed (geographical distribution)

Total article views: 723 (including HTML, PDF, and XML) Thereof 723 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 14 Sep 2025

Short summary

Intrusion of a hot basaltic magma into produces melting of host rocks. We simulate this process by means of 2D thermo-hydrodynamical model with account for variable magma properties. The melting front rapidly propagates into intact rocks, leading to a formation of a sharp temperature gradient on magma-rock interface. Our modelling proposes that a basaltic sill was injected during Krafla Fires, leading to formation of ~30 m of molten felsic magma that was penetrated by IDDP-1 drilling.


Total:	0
HTML:	0
PDF:	0
XML:	0