Structural dynamics of basaltic melt at mantle conditions with implications for magma oceans and superplumes

Transport properties like diffusivity and viscosity of melts dictated the evolution of the Earth’s early magma oceans. We report the structure, density, diffusivity, electrical conductivity and viscosity of a model basaltic (Ca11Mg7Al8Si22O74) melt from first-principles molecular dynamics calculations at temperatures of 2200 K (0 to 82 GPa) and 3000 K (40–70 GPa). A key finding is that, although the density and coordination numbers around Si and Al increase with pressure, the Si–O and Al–O bonds become more ionic and weaker. The temporal atomic interactions at high pressure are fluxional and fragile, making the atoms more mobile and reversing the trend in transport properties at pressures near 50 GPa. The reversed melt viscosity under lower mantle conditions allows new constraints on the timescales of the early Earth’s magma oceans and also provides the first tantalizing explanation for the horizontal deflections of superplumes at ~1000 km below the Earth’s surface.

This paper reports the properties of a molten basalt analogue up to 80 GPa by means of theoretical calculations. The major finding is that melt viscosity remains low up to the highest P investigated (80 GPa) with a non-monotonous P-evolution. Viscosity slightly decreases up to 20 GPa, then increases only to reach a maximum near 50 GPa (100 mPa.s) and then slowly decreases. These results differ from those obtained previously from theoretical calculations (Karki 2018 and refs therein) that did not report a slow decrease above 50 GPa, with viscosity values about 2-3 times higher than found here. Although experimental data on viscosity only exist up to modest P, experimental values are also 3-4 times higher than found here (i.e. below 10 GPa, cf highest T points in Sakamaki 2013). The decrease of viscosity above 50 GPa is attributed to the Si-O and Al-O bonds become ionic, explaining weakening of the melt structure and reduced viscosity. The most straightforward implication of such low viscosity is that the magma ocean must have cooled faster than previously thought. There are no molten basalts in the present day deep mantle except in ULVZs. Hence implications for the reported viscosity decrease in the solid lower mantle are quite stretched. There are a number of mistakes in the literature and geological background (see below), however that can be corrected. I am more concerned by the mismatch between the present results on diffusivity and viscosity and previous theoretical and experimental results, while structural and density data are in general agreement with previous works. Such mismatch, especially at modest P at which experiments can be done and are confirmed by other theoretical works, cast doubts on the new finding i.e. a decrease of viscosity at lower mantle conditions. For these reasons, I cannot recommend the paper for publication.
Geological background: I do not think that it is appropriate to state climate change as being governed by the properties of magmas at depth. This induces the idea of a control from the deep Earth on climate change, which is inaccurate and dangerous for a scientist to state. Super-eruptions (e.g. traps) did influence climate in the past, but there is no such effect since C O2 increase in the atmosphere from humanrelated activities. Most implications related to the presence of basalts at depth are fine however at places the manuscript implies that there is partial melt in the lower mantle which is not the case except in ULVZs where T are high enough to induce partial melting. Mantle plumes may rise from low-shear wave velocity provinces (that do have some ULVZs at their very base) but they are not partially molten until they almost reach the Earth's surface. Basalts plaid no role at all in core formation. Timing of magma ocean cristallisation has been revised since Abe 1997 down to a few Myr only (Lebrun 2013, Hamano 2013, the latter using 0.1 Pa.s for the viscosity of ultra-basic magma ocean.
Indeed, as pointed out in the paper, C a is larger and modifies melt's properties significantly such as increasing its viscosity. But more importantly, the authors should explain why their results obtained on basalt may be transferred to a magma ocean. That is possible for equations of state with minor corrections, but is more difficult for viscosity which strongly depends on the SiO2 content.
Melting curve of basalt: At 2200 K (the lowest of the 2 isotherms investigated here) basalt solidifies at 14 GPa (Hirose 1999). C ombination of experimental and theoretical error bars could explain that basalt is still molten in the present work at 2200 K however that can only hold up to 20-22 GPa, after which the melting curve steepens. This could explain the decreased diffusivities (and consequent viscosity increase) reported here. Even at 3000 K, the highest isotherm, basalt solidifies at 50 GPa (Hirose 1999, Gopal 2014) however the authors do not report viscosity below 40 GPa along the 3000 K isotherm. My question is therefore: may these calculations still have a meaning if done below the melting curve?
Reviewer #2 (Remarks to the Author): The authors performed ab initio molecular dynamics simulations to investigate the structure and ionic dynamics of a basaltic melt model, C a11Mg7Al8Si22O74, in a pressure range of 0-82GPa at 2,200K and 3,000K. The authors claimed the ionic and weak Si-O and Al-O bonds at high pressure leads to an anomalous trend change in transport properties near 50 GPa, which correlates well with the viscosity decrease in the lower mantle. This is an interesting observation but obviously the mantle is solid and does not have basaltic composition. Besides, the mantle has iron undergoing a spin state change, which the authors cite (not properly though. See below). The support for the melt origin of ULVZ is more relevant. Obviously, the variation of viscosity with melt composition would be most important in this case, especially the presence of iron and hydrogen (also carbon), which the authors mention and acknowledge should be investigated. The importance of these melt properties to understanding the properties of the asthenosphere is also relevant, as discussed on page 11. A comparison of this basaltic melt viscosity with the mantle viscosity structure is a bit of a stretch. Nevertheless, the correlation seems to be real and it might simply be related to the behavior of Si-O and Al-O bond strengths and coordination through the structure. It is, therefore, a relevant observation. There also exist some technical issues in the paper that could be improved. I would suggest the authors add more details to clarify these issues.
6) The description on Page 8 line 175-177 is confusing. Only the 3,000K data shows an anomalous maximum electrical conductivity at 50 GPa.
7) The Green-Kubo relation in Line 321 on Page 14 is not OK. What happens with the αβ indices?
8)The authors should provide the plot of <P(t)P(0)> to ensure the readers that it converges to 0 for viscosity calculations. It's usually very hard to converge pressure in the ab initio MD timescales. 9) How is the confidence interval of the viscosity data determined in Fig. 4? 10) At 2,200K, the change of viscosity is different from the change of diffusivity, while they are quite consistent at 3,000K. Is the Stock-Einstein relation valid under these conditions?
11) The authors cite papers of papers that mention the effect of the iron spin state change on the viscosity in the lower mantle. The first paper that made a case for thsi phenomenon is Wentzcovitch, R. M. et al. Anomalous compressibility of ferropericlase throughout the iron spin crossover. Proceedings of the National Academy of Sciences of the United States of America 106, 8447-8452, doi:10.1073America 106, 8447-8452, doi:10. /pnas.0812150106 (2009. Please make sure proper credit is given. 12) I am sure the authors can address most of these issues in their reply without difficulty.
Reviewer #3 (Remarks to the Author): This manuscript presents results from molecular dynamic calculations for a model basalt composition up to 80 GPa at 2200 and 3000K. The study finds that C N and viscosity generally increase with pressure, while diffusivity and electrical conductivity generally decrease with pressure at 2200K. Geophysical implications for the ~1000km viscosity anomaly, the ULVZs, and the magma ocean are briefly addressed.
In general, I think this paper is well written and presented. I am not an expert on MD calculations and so I am going to move forward with my review assuming that the calculations have been performed without issue. The main issue with the ms in its present form is that there are some critical references and data for 3000K pressure points are missing. These should be relatively minor edits to address. Also, the discussion of the ULVZ should be removed and the discussion for the other geophysical implications should be expanded, as this will make the paper of broader interest for the readership of Nature Geoscience. C omments: • Figure 1: thin lines should be added to guide the eye for the C N trends, similar to what is shown in Fig. S6. • Ll.101-102: The data observation or reference to the more ionic nature should be clarified. The previous discussion is about bond lengths, C N and angles, but the connection to the changing nature of the bonding is unclear, especially to the general readership of Nature Geosciences. • Ll.129, 164, 170, and others: The references need to be more diligently attended to. There are several sentences in this manuscript that talk about predictions, expectations, etc. but there is no corresponding reference. This needs to be remedied before publication. • Ll. 137 -The data below 38 GPa should be included in the supplemental material. • Ll.157 -the low pressure data for the 3000K are missing? Were the calculations at 3000K only performed from 40-70GPa? This is unclear to me in the reading of the manuscript, perhaps I just missed it, but in the methods section (ll. 300) it is implied that there is data from 0-82 GPa at 3000K. • Ll. 161 -also the diffusion coefficient of oxygen and Si increases in figure 3a (and b) in the 50-70 GPa region. • Ll. 162 -what do the authors mean that the increase is not an artifact of the calculations? The reasoning for this should be discussed at least in the supplemental data. • Figure 4: Why do diffusivity/EC and viscosity show opposite trends from 20-50 GPa? Why does diffusivity/EC data at 2200 and 3000K show different trends from 40-70 GPa, while viscosity does not? Also, what are the 'R' in the legend of figure 4? It might be good to clarify this in the figure caption.
• Ll. 263 -the discussion of the ULVZ is beyond scope of the work of this paper (~50 GPa higher pressure and 1000K hotter). This paragraph should be removed and the discussion of the 1000km viscosity anomaly, the implications for magma ocean dynamics, and/or the interpretation of geophysical observables at the mantle transition zone should be expanded.

Reviewer #1 (Remarks to the Author):
This paper reports the properties of a molten basalt analogue up to 80 GPa by means of theoretical calculations. The major finding is that melt viscosity remains low up to the highest P investigated (80 GPa) with a non-monotonous P-evolution. Viscosity slightly decreases up to 20 GPa, then increases only to reach a maximum near 50 GPa (100 mPa.s) and then slowly decreases. These results differ from those obtained previously from theoretical calculations (Karki 2018 and refs therein) that did not report a slow decrease above 50 GPa, with viscosity values about 2-3 times higher than found here. Although experimental data on viscosity only exist up to modest P, experimental values are also 3-4 times higher than found here (i.e. below 10 GPa, cf highest T points in Sakamaki 2013). The decrease of viscosity above 50 GPa is attributed to the Si-O and Al-O bonds become ionic, explaining weakening of the melt structure and reduced viscosity. The most straightforward implication of such low viscosity is that the magma ocean must have cooled faster than previously thought. There are no molten basalts in the present day deep mantle except in ULVZs. Hence implications for the reported viscosity decrease in the solid lower mantle are quite stretched. There are a number of mistakes in the literature and geological background (see below), however that can be corrected. I am more concerned by the mismatch between the present results on diffusivity and viscosity and previous theoretical and experimental results, while structural and density data are in general agreement with previous works. Such mismatch, especially at modest P at which experiments can be done and are confirmed by other theoretical works, cast doubts on the new finding i.e. a decrease of viscosity at lower mantle conditions. For these reasons, I cannot recommend the paper for publication. Geological background: I do not think that it is appropriate to state climate change as being governed by the properties of magmas at depth. This induces the idea of a control from the deep Earth on climate change, which is inaccurate and dangerous for a scientist to state. Super-eruptions (e.g. traps) did influence climate in the past, but there is no such effect since CO 2 increase in the atmosphere from human-related activities.
Reply: The calculation of transport properties requires a large model system and very long simulations. Even in the recent past, ab initio calculations could only be performed on a small system (usually < 100 atoms). In this study we expand significantly the number of atoms and carefully monitor the convergence of the calculated properties. An example on the convergence of the calculated viscosity is shown in Fig. S12 in the supplementary information. The First principles results reported here cannot be compared with calculations performed with empirical potentials particularly at very high pressures. Moreover, the trajectories from older ab initio calculations were often shorter, within 10's of ps, and the results may be not reliable. These factors may be potential reasons for the discrepancies. The reviewer has raised concern about discrepancy between our results with that of experiments. This statement is very vague as it does not mention which theoretical and experimental results. As it is that not much data on the diffusivity and viscosity of basalt exists in the literature, therefore, for decades researchers have resorted to approximate basalt using different models, starting from silicate melts, aluminosilicate melts, model basalt, etc. Thus, it is obvious that there will be differences. However, we would like to point out that this misleading remark by the reviewer is not completely fair to our work as it can severely skew the editorial decision. In fact, our results are in quite close agreement and in the same order of magnitude as that of related works, both theoretical and experimental. Some representative literature that we are obliged to cite here are the following. The viscosities reported in these works at the similar pressures and temperatures are in quite good agreement with our data.  298-315 (2018).
We agree with the reviewer and have removed the sentence insinuating climate change being a manifestation of the properties of magma at depths.
Most implications related to the presence of basalts at depth are fine however at places the manuscript implies that there is partial melt in the lower mantle which is not the case except in ULVZs where T are high enough to induce partial melting. Mantle plumes may rise from lowshear wave velocity provinces (that do have some ULVZs at their very base) but they are not partially molten until they almost reach the Earth's surface. Basalts played no role at all in core formation. Other comments:

Reply
Choice of composition: the authors explain that the very high Ca abundance of their composition is the consequence of the simplification (no Fe). Usually, petrologists compensate Fe for Mg, not Ca. Indeed, as pointed out in the paper, Ca is larger and modifies melt's properties significantly such as increasing its viscosity. But more importantly, the authors should explain why their results obtained on basalt may be transferred to a magma ocean. That is possible for equations of state with minor corrections, but is more difficult for viscosity which strongly depends on the SiO2 content. Melting curve of basalt: At 2200 K (the lowest of the 2 isotherms investigated here) basalt solidifies at 14 GPa (Hirose 1999). Combination of experimental and theoretical error bars could explain that basalt is still molten in the present work at 2200 K however that can only hold up to 20-22 GPa, after which the melting curve steepens. This could explain the decreased diffusivities (and consequent viscosity increase) reported here. Even at 3000 K, the highest isotherm, basalt solidifies at 50 GPa (Hirose 1999, Gopal 2014) however the authors do not report viscosity below 40 GPa along the 3000 K isotherm. My question is therefore: may these calculations still have a meaning if done below the melting curve?  (2007)].

Reviewer #2 (Remarks to the Author):
The authors performed ab initio molecular dynamics simulations to investigate the structure and ionic dynamics of a basaltic melt model, Ca 11 Mg 7 Al 8 Si 22 O 74 , in a pressure range of 0-82GPa at 2,200K and 3,000K. The authors claimed the ionic and weak Si-O and Al-O bonds at high pressure leads to an anomalous trend change in transport properties near 50 GPa, which correlates well with the viscosity decrease in the lower mantle. This is an interesting observation but obviously the mantle is solid and does not have basaltic composition. Besides, the mantle has iron undergoing a spin state change, which the authors cite (not properly though. See below). The support for the melt origin of ULVZ is more relevant. Obviously, the variation of viscosity with melt composition would be most important in this case, especially the presence of iron and hydrogen (also carbon), which the authors mention and acknowledge should be investigated. The importance of these melt properties to understanding the properties of the asthenosphere is also relevant, as discussed on page 11. A comparison of this basaltic melt viscosity with the mantle viscosity structure is a bit of a stretch. Nevertheless, the correlation seems to be real and it might simply be related to the behavior of Si-O and Al-O bond strengths and coordination through the structure. It is, therefore, a relevant observation. There also exist some technical issues in the paper that could be improved. I would suggest the authors add more details to clarify these issues.
1) The authors should clarify the definition of CN in the current liquids. It can be based on a sharp bond-length threshold, or integration of the first peak in g(r), or Voronoi tessellation.

Reply: We thank the reviewer for pointing this out and have included the description in the methods section of the manuscript also. The CN has been calculated by integration of the first nearest neighbour peak of the radial distribution function, g(r), up to the first minimum.
2) The authors described bond angle changes in Page 4-5. It would be more informative to show a typical bond angle analysis, instead of the plain description.

Reply: On request from the reviewer, we have included bond angle distribution plots for O-Si-O and O-Al-O for different pressure points at 2200 K, in the supplementary section. The ~90 o is the octahedral O-Si-O and the ~170 o is the axial (linear) O-Si-O showing that the local environment is almost octahedral. The O-Al-O becomes 6 coordinated much faster than O-Si-O. This is also
reflected in the coordination number analysis. This has also been written in the manuscript now and we hope will assist in more lucid reading of the manuscript by the readers.
4) The discussion of diffusivity change is a bit confusing and might be improved. The authors claim Mg and Ca move easily because there are open spaces. But the diffusivity data in Fig.3 shows Ca or Mg are usually the slowest elements under most higher pressures at 2,200K.
Reply: At 60 and 70 GPa, the diffusivities have been recalculated and definitely more plausible magnitudes have been reported now. We apologize for the error.

5) At 60 GPa and 70
GPa, Mg becomes the slowest specie while it is fastest at 50 GPa and at 80 GPa. Why does the trend change so rapidly? I suggest the authors include the confidence interval for the diffusivity data.
Reply: At 60 and 70 GPa, the diffusivities have been recalculated and definitely more plausible magnitudes have been reported now. We apologize for the error. It is very difficult to show confidence levels in log plots. So we have written the confidence levels in the manuscript.
6) The description on Page 8 line 175-177 is confusing. Only the 3,000K data shows an anomalous maximum electrical conductivity at 50 GPa.
Reply: We have added a line in the relevant section that for 2200 K, an anomalous maximum is also seen at 60 GPa. Hence, anomalous maxima are seen for both 2200 and 3000 K. This should remove any confusion.
7) The Green-Kubo relation in Line 321 on Page 14 is not OK. What happens with the αβ indices?
Reply: The Green-Kubo equation had a factor of 3 missing in the denominator. It was a typo and has now been corrected. We thank the referee for pointing this out. The meaning of αβ have also been explained in the revision. The indices refer to the three off-diagonal elements of the stress tensor, namely, xy. yz, and zx components. Therefore, the theoretical background and reproducibility of our calculations has been ensured.
8) The authors should provide the plot of to ensure the readers that it converges to 0 for viscosity calculations. It's usually very hard to converge pressure in the ab initio MD timescales.
Reply: Within the time scale of ab initio MD, sometimes it is indeed difficult to converge the stress autocorrelation function (SACF) to 0. Hence, a plot showing the decay of the SACF and convergence to zero has been added in the supplementary section (Fig. S11). This valuable suggestion definitely makes the viscosity analysis more convincing.
9) How is the confidence interval of the viscosity data determined in Fig. 4?
Reply: We divided the MD trajectory into several time segments with different time origins to compute the correlation functions while making sure that convergence had been achieved in each segment as described above. The error was the standard deviation from the mean value. We have now explained this in the methods section.
10) At 2,200K, the change of viscosity is different from the change of diffusivity, while they are quite consistent at 3,000K. Is the Stock-Einstein relation valid under these conditions?
Reply: The Stokes-Einstein relation between particle diffusivity and fluid viscosity works well for comparatively simpler liquids such as metals [Poirier, J. P., Transport properties of liquid metals and viscosity of the Earth's core. Geophys. J. 92, 99-105, (1988);Dobson, D. P., et al., In situ measurement of viscosity of liquids in the Fe-FeS system at high pressures and temperatures. Am Mineral 85, 1838-1842(2000]. However, there have been reports about discrepancies for fluids of more complex composition implying that the diffusion involves mechanisms that cannot be explained by such a simple model [Rev. Mineral. Geochem. 72, 971-996 (2010)]. In the Stokes-Einstein relation, the effective hydrodynamic radius of the atoms has to be either taken as an average of the different species involved or that of the atoms which contribute majorly to the diffusivity. This can give rise to significant errors and therefore the Green-Kubo method of finding out the viscosity from MD is definitely more reliable.
11) The authors cite papers of papers that mention the effect of the iron spin state change on the viscosity in the lower mantle. Reply: This is a justified claim. In the revised manuscript, as requested by this reviewer, we have removed the major discussion on the ULVZ and the geological implications of our results on it.

Reviewer #3 (Remarks to the Author):
This manuscript presents results from molecular dynamic calculations for a model basalt composition up to 80 GPa at 2200 and 3000K. The study finds that CN and viscosity generally increase with pressure, while diffusivity and electrical conductivity generally decrease with pressure at 2200K. Geophysical implications for the ~1000km viscosity anomaly, the ULVZs, and the magma ocean are briefly addressed.
In general, I think this paper is well written and presented. I am not an expert on MD calculations and so I am going to move forward with my review assuming that the calculations have been performed without issue. The main issue with the ms in its present form is that there are some critical references and data for 3000K pressure points are missing. These should be relatively minor edits to address. Also, the discussion of the ULVZ should be removed and the discussion for the other geophysical implications should be expanded, as this will make the paper of broader interest for the readership of Nature Geoscience. Comments: • Figure 1: thin lines should be added to guide the eye for the CN trends, similar to what is shown in Fig. S6. Fig. 1 for better aid to the eye.

Reply: We have now joined the points in
• Ll.101-102: The data observation or reference to the more ionic nature should be clarified. The previous discussion is about bond lengths, CN and angles, but the connection to the changing nature of the bonding is unclear, especially to the general readership of Nature Geosciences.
Reply: When the coordination changes from four fold to six fold, bond lengths increase to accommodate more atoms in the polyhedron. Therefore, the bonds become weaker going through an eventual transition from covalent (strong) to ionic (weak). We have added a line on this in the manuscript too.
• Ll.129, 164, 170, and others: The references need to be more diligently attended to. There are several sentences in this manuscript that talk about predictions, expectations, etc. but there is no corresponding reference. This needs to be remedied before publication.
Reply: We agree with the reviewer and have now cited the appropriate references accordingly.
• Ll. 137 -The data below 38 GPa should be included in the supplemental material.
Reply: We have kept this for better representation to show the low pressure tendencies as well. This is important to show the change in the structural properties that lead to the change in the transport properties.
• Ll.157 -the low pressure data for the 3000K are missing? Were the calculations at 3000K only performed from 40-70GPa? This is unclear to me in the reading of the manuscript, perhaps I just missed it, but in the methods section (ll. 300) it is implied that there is data from 0-82 GPa at 3000K. • Ll. 263 -the discussion of the ULVZ is beyond scope of the work of this paper (~50 GPa higher pressure and 1000K hotter). This paragraph should be removed and the discussion of the 1000km viscosity anomaly, the implications for magma ocean dynamics, and/or the interpretation of geophysical observables at the mantle transition zone should be expanded.

Reply:
We have now added the results and a discussion on an addition calculation performed at at 120 GPa and 4000 K which is directly relevant to the discussion of the ULVZ. Also, we have highlighted the viscosity anomaly at ~50 GPa providing the first tantalizing explanation for the horizontal deflections of superplumes from the core-mantle boundary.
The revised manuscript addressed all my comments satisfactorily. I also see that the authors went a long way to address the concerns of other referees. From my perspective, this is a well presented and thoughtful piece of research reporting a curious drop in viscosity of a basaltic melt above 50 GPa. Geophysical implications are not always obvious and can be argued about, but the results seem credible, and that is more important to me. The latter will stand, while the former may not, depending on progress in understanding in the field. I believe the paper is suitable for publication.
Reviewer #3 (Remarks to the Author): This study investigates the effect of pressure and temperature on the transport properties of basaltic melt model composition C a11Mg7Al8Si22O74, at pressures from 0-120 GPa and in the temperature range of 2200-4000K (not all pressures and temperatures combined).
The added supplemental figures (i.e. S1b, S1d, S2b, S1d) have improved the comparison to previous data and show many of the same general trends, which was a concern of all the previous reviewers. The authors also found a calculation error that has been corrected (e.g. Reviewer 2 point 4). Also, the focus of the discussion has been shifted and now focuses on plumes and the ULVZ. The applicability of the data to Generally, I feel that the requested issue that were brought up by the previous reviewers have been addressed satisfactorily. I have two outstanding concerns prior to publication. They are given below. There are a few minor typos (e.g. ll. 193), I would consider a final edit to clean up the ms.
C omments: 1) Focusing on the ULVZ and not having Fe present in the melt does create some questions as the ULVZs are thought to possibly be Fe rich -they are sitting atop the core -a giant vat of Fe (and other elements). More over while the 120 GPa data point is at significantly higher pressure than the previous versions data -which had a maximum pressure of 82 GPa -the C MB is at 135 GPa and the ULVZs are only 10's on km high -meaning that the data calculated here for 120 GPa is still 100's of kms (approaching 500km) too shallow. C onsidering these two issues (lack of Fe and low pressure), I would consider still removing this discussion and focusing on the plume story. In fact, I would consider removing the 4000K, 120 GPa data point all together and save it for a followup study or something. As shown in Figure 4, there is every little that can be said about viscosity at the C MB as there is no trend present (i.e. only one point). There could be a follow up section in the discussion like "If this trend continues with depth and temperature, there may be implications for the ULVZ by . . . ." but anything more is really too speculative with the data present.
2) The added figures in this version of the manuscript are nice. I would suggest adding the diffusivity and viscosity data of others (both calculations and experiments) to Figure 3 and 4. I.e. refs 19-21, 2-4, etc. At a minimum this type of plot or table should be included in the supplemental material. As long as this data is generally consistent, I would think that this manuscript is ready for publication.