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

Gerbault, Muriel; Melnik, Oleg; Borisova, Anastassia

doi:10.5194/se-16-1307-2025

Articles | Volume 16, issue 11

https://doi.org/10.5194/se-16-1307-2025

© Author(s) 2025. This work is distributed under
the Creative Commons Attribution 4.0 License.

https://doi.org/10.5194/se-16-1307-2025

© Author(s) 2025. This work is distributed under
the Creative Commons Attribution 4.0 License.

Articles | Volume 16, issue 11

Research article

|

12 Nov 2025

Research article |

| 12 Nov 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

Download

Final revised paper (published on 12 Nov 2025)
Preprint (discussion started on 07 Jul 2025)

Interactive discussion

Status: closed

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
- AC3: 'Reply on EC1', Oleg Melnik, 21 Oct 2025
  
  We have submitted the corrected version of the manuscript at the end of September.
  
  Citation: https://doi.org/10.5194/egusphere-2025-3036-AC3

Peer review completion

AR: Author's response | RR: Referee report | ED: Editor decision | EF: Editorial file upload

AR by Oleg Melnik on behalf of the Authors (25 Sep 2025) Author's response Author's tracked changes Manuscript

ED: Publish subject to minor revisions (review by editor) (30 Sep 2025) by Virginie Pinel

AR by Oleg Melnik on behalf of the Authors (01 Oct 2025) Author's response Author's tracked changes Manuscript

ED: Publish as is (02 Oct 2025) by Virginie Pinel

ED: Publish as is (21 Oct 2025) by Andrea Di Muro (Executive editor)

AR by Oleg Melnik on behalf of the Authors (27 Oct 2025) Manuscript

Short summary

We present a 2D numerical modelling of magma-rock interaction after a sill injection into shallow crust. A solver based on OpenFOAM package simulates heat and momentum transfer between magmas of different compositions, with contrasting properties. Two scenarios are considered to explain the presence of fresh rhyolitic fragments excavated during drilling at Krafla caldera in 2009: partial melting of felsic crust triggered by either a 300 m thick rhyolite intrusion or a 100 m thick basalt sill.