硅酸鹽熔體和液態鐵的兩相模擬與核幔間的氧交換

張志剛， 劉錦， 張毅剛， 李娟

1 中國科學院地質與地球物理研究所， 中國科學院地球與行星物理重點實驗室， 北京 100029 2 中國科學院大學地球與行星科學學院， 北京 100049 3 北京高壓科學研究中心， 北京 100094 4 中國科學院深地科學卓越創新中心， 廣州 510640

0 引言

和太陽系其他行星一樣，地球起源于星際塵埃，通過不斷吸積、碰撞，經歷星子、星胚等階段，在大約一億年時間里逐漸長大成原始地球，而后又很可能遭到火星大小的星球引發的大撞擊事件，形成了地月系統，此后逐漸冷卻至今(Stevenson, 2008).在此過程中，地幔中許多元素隨地球生長而被帶入地核，而后又隨著地球冷卻被緩慢釋放回地幔.核幔之間的這種元素交換過程伴隨著超長時空尺度的物質和能量流動，是地球成為宜居星球的重要原因(Langmuir and Broecker, 2012).

氧是地幔中豐度最高的元素之一，也是地核中重要的化學組分，因而是解譯地球演化歷史的關鍵線索.實際上，氧元素廣泛參與了核幔間其他諸多元素的交換，與這些元素在地幔和地核中的分布有很強的相關性(Chidester et al., 2022).氧在核幔間的交換與地幔中的FeO含量和氧化還原狀態緊密聯系，因而是地幔中巖石和地球化學過程的主控因素之一(Wade and Wood, 2005).此外，最近華人科學家提出和主導了“地幔底部富氧假說”(Mao and Mao, 2020)，讓人們把下地幔底部的地震學觀測異常(Garnero et al., 2016; He et al., 2021)與富氧礦物聯系起來(Liu et al., 2017)，這些富氧礦物與地核間的氧交換也是未來值得研究的重要問題.

氧元素在核幔間交換的底層機制是其在核幔分異過程中的平衡(Wood et al., 2006).為定量刻畫這些平衡，在過去的二三十年時間里人們開展了大量的實驗和建模工作(Chidester et al., 2022).但已有研究仍然存在數據溫壓區間有限、模型差別大的問題.本研究我們運用第一性原理計算，構建了能夠準確預測硅酸鹽熔體和液態鐵之間元素交換的方法，基于模擬結果和已有高溫高壓實驗數據建立了新的平衡模型，最后運用這些結果探討了核幔間的氧交換與地核的氧含量等問題.

1 研究方法

1.1 模擬溫壓條件的選擇

如圖1所示，我們將前人所做的近兩百個核幔平衡實驗點投到溫壓分布圖上(具體數據值和引文詳見Chidester等(2022)的附錄)，可以看到大部分數據點集中在4000 K、80 GPa以下，目前最高溫壓的實驗是Fischer等(2015)提供的，達到5700 K、100 GPa.

圖1 硅酸鹽熔體與鐵合金液體元素分配實驗的溫壓范圍 圖中藍線和紅線分別為橄欖巖成分地幔的 固相線和液相線(Fiquet et al., 2010).Fig.1 T-P distributions of the existing experimental data on element partitioning The blue and red curves are the experimentally determined solidus and liquidus respectively for the mantle with a peridotite composition (Fiquet et al., 2010).

圖1中還給出了橄欖巖成分地幔的固相線和液相線(Fiquet et al., 2010)，顯然已有實驗數據的溫度基本位于固相線以上，并且很大一部分落在地幔固液兩相區間.這主要是由于現代巖漿洋理論認為：地球早期大型撞擊熔融所致的金屬鐵液滴，與全球尺度的硅酸鹽熔體(即所謂巖漿洋)充分混合，在巖漿洋的底部達到化學平衡；如果巖漿洋與原始地核之間有固結的地幔，這些金屬鐵液滴將在重力作用下快速穿過固體地幔與原始地核融合，所以地核所能獲取的氧和其他元素應當主要決定于巖漿洋底部的平衡(Wade and Wood, 2005).

最近幾年的研究表明，地球早期巖漿洋的深度和溫度也許要遠超預期.比如，地月系統可能形成于更大力度的星球撞擊事件，大撞擊所釋放出的巨大能量足以將整個地幔全部熔融(Lock et al., 2020).如此規模的巖漿洋在冷卻過程中，也許會首先從地幔中部固結，從而有利于地幔深處更長時間地保留硅酸鹽熔體，使得核幔邊界處氧的交換得以更有效地進行(Davies et al., 2020).

為了能夠探討與超深巖漿洋假說相關的核幔氧交換作用，本研究主要考察了135 GPa條件下的4000 K、5000 K和6000 K三個溫度點.如圖1所示，這些溫壓條件的模擬是對已有實驗的有益補充.

1.2 模擬及分析方法

本研究我們結合基于密度泛函理論的第一性原理方法與機器學習方法，構建了超大體系的兩相共存模擬.根據我們前期的技術積累(Zhang et al., 2022, 2020)，這些模擬結果準確與否的關鍵在于原子間相互作用勢能的計算.與傳統的勢能模型不同，新一代的機器學習模型不再受限于固定的函數形式，能在高維度上實現對第一性原理模擬結果的準確重現和預測，與此同時模擬的效率提高至少三個數量級以上.我們在Zhang等(2022)的附錄中對此方面相關技術的細節進行了詳細闡述，這里不再贅述.

在建立了高質量的機器學習勢能模型以后，我們得以進行超過一千個原子的分子動力學模擬.在測試了2048、4096、8192個原子體系的有限尺寸效應之后，我們在模擬體系中采用了4096個原子.根據McDonough(2014)提出的地球總體成分模型，模擬盒子中設置了1968個氧原子、704個鐵原子、608個硅原子和816個鎂原子.為了得到相對規則的兩相分布以便于后續分析，借用前人在流體體系的研究經驗(Muller et al., 2020)，我們將初始鐵原子集中于模擬盒子中間(代表地核)，其余硅、氧、鎂組分都放置于與之共存的硅酸鹽熔體相中(代表地幔巖漿洋)，在經過大約100 ps(100000步)的分子動力學模擬后形成如圖2a所示的構型，此時監測的能量和結構隨模擬時長的波動顯示兩相已經過充分互溶達到平衡，而后模擬又運行了200 ps(200000步)用以統計各熱力學參量.為了測試初始構型對結果的影響，我們還另外進行了幾組模擬，在保持總體成分一定的情況下改變了兩相的初始成分，模擬結果表明，平衡后的統計結果在誤差范圍內不變.

超大體系的設定使得硅酸鹽熔體相和液態鐵相的界面較為穩定和平滑，非常有利于對模擬結果的分析.如圖2b所示，沿著與兩相接觸面垂直的Z軸方向，通過統計局部密度分布，可以很好地區分出硅酸鹽相、液態鐵相以及介于兩者之間的界面區.我們采用下面的函數對密度分布進行擬合(Muller et al., 2020)：

(1)

該式中ρsil和ρmet分別為硅酸鹽相和金屬鐵相的共存密度，z1和z2為圖中左右兩個界面的位置，δ1和δ2則為兩個界面區的厚度.

圖2 5000 K、135 GPa條件下核幔平衡的兩相共存模擬(a) 分子動力學模擬的原子構型，由4096個原子組成，其中紅色為氧、藍色為鐵、黃色為硅、粉色為鎂； (b) 沿兩相共存的Z方向密度 分布統計圖.Fig.2 Coexistence simulation for the core-mantle chemical equilibrium at 5000 K and 135 GPa(a) A typical atomic configuration of the molecular dynamics trajectory with 4096 atoms (red: oxygen; blue: iron; yellow: silicon; pink: magnesium)； (b) Density profile along the elongated Z direction.

我們進一步對Z軸方向上各元素的局部摩爾分數進行分析.如圖3a所示，這些分布圖與密度分布類似，可以根據由密度分布所劃分的相區分別進行原子比例統計，這使得我們可以排除界面區表面張力的影響，從而更為準確地得到兩相共存時的成分.利用統計得到的兩相密度和成分，我們還分別對硅酸鹽相和液態鐵相進行了單相分子動力學模擬，得到了徑向分布函數和配位數等微觀結構信息，如圖3b所示.

1.3 平衡模型

前人的實驗和我們開展的計算模擬，為核幔間氧和其他元素的交換提供了約束，但其涉及的溫壓和成分總體還很有限.在實驗和模擬數據的基礎上發展相關的熱力學模型，不僅能更深入地理解核幔間的平衡，對于核幔演化的相關討論也是非常必要的.

圖3 5000 K、135 GPa條件下核幔平衡模擬中，(a) 沿兩相共存的Z方向各元素摩爾分數(x)統計圖； (b) 氧元素與其他元素配對的徑向(r)分布函數g(r)結構分析Fig.3 (a) Mole fractions (x) of four types of atoms along the Z direction and (b) radial (r) distribution functions g(r) related with oxygen for the coexisting simulation under 5000 K, 135 GPa

地幔中的氧通常以陰離子形式與不同的陽離子結合成硅酸鹽或者氧化物，而地核中的氧則一般與金屬鐵鎳形成零價的合金溶液，因此氧在核幔間的交換實際上是與陽離子的交換分不開的(Wade and Wood, 2005).對于價態為+n的陽離子M，由于地幔中氧的價態一般為-2，氧伴隨M進出地核的平衡可以表達成下面的“溶解”平衡(這里稱之為(R1))：

式中M和O表示地核金屬鐵合金中的組分，MOn/2表示地幔中的氧化物組分.

如果將FeO的溶解平衡式與其他MOn/2氧化物的溶解平衡式聯立，可以得到以下“交換”平衡(這里稱之為(R2))：

從平衡原理上看，(R1)和(R2)兩種表征方法是等價的，但實際應用時，基于各種假設和簡化，不同平衡式對實驗數據回歸的總體符合程度會有所差異(Badro et al., 2018; Chidester et al., 2022).我們在本研究中僅采用溶解平衡(R1).

(2)

該式中的下標met和sil分別表示金屬鐵相和硅酸鹽相.

(3)

式中T、P為溫度和壓力，ΔG0為該溫度壓力下純氧化物和金屬的自由能差，可以進一步在第二個等號右邊近似表達為標準態的焓變ΔH0、熵變ΔS0和體積變化ΔV0的關系.而根據前人研究，ΔH0、ΔS0和ΔV0一般隨溫壓變化較小，因此可以將與之對應的a,b,c近似地設為常數(Wade and Wood, 2005).

-lnγM,met，

(4)

式中γ為液態鐵和硅酸鹽熔體中各組分的活度系數.

在實際建立式(4)的模型時，可以進一步簡化以減少待解參數的數目.比如，硅酸鹽熔體中各氧化物的活度系數γMOn/2,sil嚴格來說應該是隨溫壓和組成變化而變化的(Mukhopadhyay et al., 1993)，但前人研究發現這一項對結果的影響總體較小；另外，液態鐵相中鐵的活度系數γFe,met總體的變化范圍也不大(Fischer et al., 2015; Helffrich et al., 2020; Liu et al., 2020).所以，可將式(4)中這兩部分的影響合并到前面的參數a、b、c中去.

對于式(4)平衡模型剩下的液態鐵中各元素的活度系數γM,met和γO,met，可采用Ma (2001)提出的模型描述：如果將除鐵之外的元素記為i，其活度系數為

(5)

(6)

2 計算結果

2.1 模擬結果

我們對5000 K、135 GPa下的兩相結構進行了分析.如圖3b，結果表明四種元素從硅酸鹽相到液態鐵相的微觀結構發生了很大變化，其中Si-O、Mg-O的徑向分布函數變化最為明顯.由于兩相中的元素成分差異大(圖3a)，各元素周圍的氧配位數也發生較大變化：在硅酸鹽相中，O、Mg、Fe、Si周圍的氧配位數分別為10.5、8.0、5.7和6.0；而在液態鐵相中，這些氧配位數分別為1.9、4.2、1.1和0.7.這些結果與前人對下地幔硅酸鹽熔體和地核鐵合金微觀結構分析結果基本一致(Alfè et al., 1999; de Koker et al., 2013).

圖4為三個溫度下兩相模擬所得到的平衡成分.從中可以看出，氧、硅、鎂在地核中的含量隨溫度快速增加，6000 K時氧和硅大量進入地核，可以分別達到26.2 mol%和17.5 mol%，即使是高度親石的鎂在地核中也可以達到4.1 mol%.前人對地月系統形成過程的數值模擬表明，大撞擊導致的核幔溫度甚至可以超過7000 K(Canup, 2012; Lock et al., 2020)，因此地球早期進入地核的元素含量在特定階段很可能超過圖4a中的估計.Wahl和Militzer (2015)通過估算Fe-O-Mg體系的自由能，提出在核幔邊界壓力下、大約7000 K就足以使得地幔與地核物質完全混溶，我們的結果雖然在趨勢上一致，但要達到Fe-O-Si-Mg體系的完全混溶溫度顯然需要比7000 K高很多.

隨著核幔邊界溫度降低至5000 K，地核中氧和鎂的含量迅速降低至12.4 mol%和0.6 mol%，而硅的含量則穩步下降至12.6 mol%；而到4000 K時，僅有2.8 mol%的氧和0.1 mol%的鎂還留在地核，此時地核中主要的輕元素為硅，含量仍可達到9.0 mol%.

與地核成分隨溫度快速變化不同，圖4b顯示平衡時地幔的成分總體變化不大，FeO、SiO2和MgO端元的濃度大致為2.3 mol%、38.0 mol%和59.7 mol%左右.

我們在圖4中還示意了前人綜合地球物理觀測、地球化學分析和礦物物理計算所估計的現今地核與地幔成分(Badro et al., 2014; Workman and Hart, 2005).由于地幔完全固結會在很大程度上阻礙元素在核幔的交換平衡(Davies et al., 2018)，因此如果地球早期巖漿洋深達核幔邊界，現今地核與地幔的成分應當在核幔邊界溫度降至固相線和液相線(根據圖1中(Fiquet et al., 2010)的實驗分別為4180 K、5375 K)對應的溫度區間時保留下來.將這些成分估計分別畫在圖4a和4b中，可以看出我們模擬所采用的McDonough(2014)地球總體成分模型將形成一個貧氧的地核和一個更為還原的地幔.如果要與現今地核與地幔的成分吻合，整個地球體系中需要加入更多的氧，我們將在后文結合平衡模型做進一步討論.

為了與圖1中涉及到的實驗數據進行比較，我們在圖5中給出了用式(2)計算得到的三種氧化物端元溶解平衡的交換系數.可以看到盡管實驗數據分布較散，在2500 K以上的數據呈現出隨溫度升高而增大的趨勢.我們的模擬結果與此相符，并基本落在實驗結果的分布范圍內，MgO的交換系數落在實驗數據分布的邊緣，表明模擬可能在一定程度上低估了地核中鎂的含量.

由式(4)可知，交換系數的數值不僅與溫壓有關，還與成分有很大的關聯，這是造成圖5中交換系數分散的重要原因.利用我們下面建立的平衡模型，可以進一步考察實驗結果與模擬結果的異同.

圖4 兩相共存模擬得到的135 GPa和不同溫度下液態鐵相中的(a)O、Si、Mg摩爾分數和(b)硅酸鹽熔體相中的FeO、SiO2、MgO摩爾分數.圖中顏色與圖例相對應的幾條水平線代表現今地核與地幔的估計成分(Badro et al., 2014; Workman and Hart, 2005)Fig.4 Mole fractions of O, Si and Mg in the iron phase(a) and FeO, SiO2 and MgO mole fractions in the silicate melts (b) from the coexisting simulations under 135 GPa and three temperatures. The horizontal lines with colors noted in the legends are the estimated compositions of current core and mantle by previous studies (Badro et al., 2014; Workman and Hart, 2005)

圖5 實驗和計算模擬得到的氧化物溶解平衡交換系數Fig.5 Exchange coefficients for dissolving reactions of oxides from experiments and simulations

2.2 平衡模型

我們將Fe-O-Si-Mg體系兩相共存模擬表達成FeO、SiO2和MgO的三個獨立平衡式(如表1所列)，這三個平衡分別有3個與溫壓有關的參數，而液態鐵中輕元素的相互作用參數另有6個獨立參數，因此總共有15個待解參數.

圖1和圖5中的實驗數據點大多數都有非常復雜的組分，其中不少實驗的硫、碳含量甚至超過40%，這些組分可能對本研究關注的氧、硅、鎂的交換有比較重要的影響(Fischer et al., 2015).另外，根據前人的研究，核幔分異時決定核幔成分的巖漿洋底部壓力至少在35 GPa以上(Badro et al., 2015; Fischer et al., 2015).因此，我們在建模時只采用30 GPa以上、并且硫、碳含量不超過5%的數據點用于參數擬合.這樣篩選后，三個平衡共有73個數據點可用于擬合.

平衡模型的參數可以通過線性最小二乘法擬合得到，此過程中我們對擬合參數進行了統計學顯著性F檢驗，在P值0.05水平進一步剔除了6個統計意義不明顯的參數，最終得到9個獨立參數的擬合值(見表1).如圖6，盡管該模型做了很多簡化，三個平衡交換系數的模型計算值與測量值總體較為符合：擬合的R2值為0.98，最大偏差約為1.9個自然對數單位(相當于0.82個以10為底的對數單位)，95%置信區間的均方根誤差為0.63個自然對數單位(相當于0.27個以10為底的對數單位)，與實驗的誤差相當(Chidester et al., 2022; Du et al., 2019;Fischer et al., 2015).

表1 Fe-O-Si-Mg核幔平衡模型的參數擬合值Table 1 Regressed parameters for the core-mantle reaction equilibria of Fe-O-Si-Mg systems

圖6 交換系數自然對數的測量值與模型預測值的對比 該圖僅采用了圖1和圖5中壓力大于30 GPa、并且硫和碳摩爾分數小于0.05的數據點.Fig.6 Natural logarithm values of exchange coefficients from measurements and model predictions Here we only used subset of data points at pressures higher than 30 GPa in Fig.1 and 5 with mole fractions of sulfur and carbon less than 0.05.

上述平衡模型能較為準確地同時重現高溫高壓實驗和本研究的計算模擬結果，在一定程度上也反映了我們構建的兩相共存模擬可以比較準確地預測核幔間的元素交換.盡管如此，我們也注意到4000 K時SiO2溶解平衡的交換系數仍有可觀的偏差.由圖1我們可以看到，該溫度略低于地幔的固相線溫度(4180 K)，因此也許需要更長時間的模擬以使得系統達到真正的平衡.

3 討論與啟示

現今地核中的氧含量有多少？又是如何通過核幔間的氧交換演化至今的？這兩個相互關聯的問題一直是地球深部領域富有爭議的重要研究方向.我們在前期的論文中對此有過初步的分析和討論(Zhang et al., 2022).借助本研究新建立的核幔平衡模型，我們將做更進一步的討論.

傳統的地球生長模型認為，地球初始成分非常還原，隨著地球的吸積、碰撞長大變得越來越氧化，最終達到現今地幔與地核的狀態(Wood et al., 2006).由于地幔的成分相對明確，人們往往通過地幔中FeO含量以及Ni、Co、Cr、V等元素在地幔中的豐度來判別生長路徑的成功與否.可想而知，由于每種生長路徑都涉及到很多不確定參數，導致人們對本節開頭提出的兩個問題有很多不同的解答.

圖7 地核中的氧和硅含量 圖中實線為本研究得到的巖漿洋底部不同壓力下與地幔巖成分共存時地核中氧和硅的含量，每條實線上的數字代表相應的壓力(單位為GPa)，沿著箭頭的方向從固相線溫度升高到液相線溫度.左邊陰影區代表地球物理觀測和礦物物理約束的地核可能的氧和硅的含量區間，陰影區中間的五角星代表地核成分的最優解.圖中其余數據點代表前人用不同地球生 長演化路徑得到的地核成分.Fig.7 Contents of oxygen and silicon in the Earth′s core The solid curves are those coexisting with the pyrolite composition silicate melt under different temperatures and pressures. The pressure of each curve (in GPa) is denoted with a number on it and the temperature is increasing from the solvus to the liquidus temperature along the curve with the direction noted by the arrow. The shadowed regime is the most probable range of oxygen and silicon in the current core as determined by the seismic observations and mineral physics, within it the star would be the best solution. The other points in theFigure are the compositions estimated by various studies with different accretion paths.

我們認為將今論古是約束上述問題的更為有效的途徑.首先，現今地核中究竟有多少氧，雖然仍存在爭議，但從地球物理觀測和礦物物理學等諸多研究來看，似乎地核應該相對富氧，并且這些氧主要富集在液態外核之中.這方面代表性的工作有：英國倫敦大學學院的Dario Alfe等通過研究氧、硅、硫在內外核的分配系數，發現氧是造成內外核密度差、驅動地球磁場的主要原因，其含量可以達到13.0 mol%(～4.3 wt%)(Alfè et al., 2002, 2000).法國巴黎地球物理學院的James Badro聯合英國倫敦大學學院的John Brodholt等，基于二元體系和理想混合規則，得到了內外核邊界和核幔邊界處(Fe-Ni)-C-O-Si-S液體合金的密度和波速，通過與地球物理觀測數據進行比對，他們發現氧是外核中不可或缺的輕元素，其含量的最優解大約為3.7 wt%，最高可達到5.8 wt%(Badro et al., 2015, 2014).最近，Umemoto 和Hirose(2020)將這種思路拓展到(Fe-Ni)-C-O-Si-S-H液體合金體系，發現地核中氧的可能含量還與地核中的氫以及地核溫度相關聯：相對高溫的地核中(內外核邊界溫度為6000 K)氧含量可以達到6.6 wt%.我們在圖7中將上述研究確定的氧和硅含量范圍用陰影區表示出來.其中的黑色五角星是Badro等(2014)所得到的地核成分“最優解”.

地核又是如何獲得相對高的氧含量呢？這同樣可以以終為始往前推論：如前所述，地核的成分決定于巖漿洋底部溫壓條件下的核幔平衡，因為如果地幔開始結晶至最后完全固結，核幔間元素交換的效率會越來越低；如果假設與地核共存的地幔成分固定為現今的地幔巖模型成分(FeO、SiO2和MgO分別大約為6.2 mol%、38.7 mol%和49.9 mol%(Workman and Hart, 2005))，我們可以利用本研究得到的模型(表1)計算出不同溫度和壓力下地核的成分.如圖7所示，我們發現地核中氧的含量會隨著壓力升高而升高，并且總體來說硅含量都要大于氧的含量，這與前人采用復雜的地核生長模擬得到的結果是一致的(Fischer et al., 2017, 2015; Rubie et al., 2015)，表明地核的成分確實主要是由最后一次巖漿洋事件決定的.

如果將圖7的實線和圈定的地核氧、硅含量區間放到一起，可以看到只有壓力大于80 GPa以上、并且相對低溫(接近相應壓力的固相線溫度)的巖漿洋底部才能形成與觀測相符的現今地核.這恰恰是前人困擾和爭議的地方，因為用地幔中Ni、Co、Cr、V的含量做判據，地核生長模型認為巖漿洋底部比較合理的壓力在40～75 GPa范圍(Badro et al., 2015; Fischer et al., 2015).

為了解決以地幔成分為判據所得到的地核成分與觀測不相符的問題，Siebert等(2013)和Badro等(2015)設計了與傳統地球生長模型不同的從氧化到還原的路徑.從圖7以及原文列出的數據我們看到，此路徑下可以在相對低壓的條件下實現地核氧的富集.其原因在于：這些研究認為巖漿洋最深只到1700 km左右，意味著地幔底部還有很大的已經固結區域；地核中的氧原本要隨著體系還原而降低含量，但受限于固結地幔使得其仍然被“鎖定”在地核之中.盡管從氧化到還原的生長路徑看似可以解決目前的矛盾，但又帶來了新的問題：地球為什么從高度氧化的初始狀態逐漸還原至今？這方面仍然沒有合理的解釋，相反的論證卻經常可以看到，比如Rubie等(2015)的系統研究就認為早期地球因為靠近太陽，應該是高度還原的.

如果地球的生長確實更應該遵循著從還原到氧化的演化路徑，也就是說地核中的氧含量隨地球生長而遞增，從圖7可知，最后一次幾乎深達核幔邊界的超深巖漿洋事件看來是不可避免的.按照前人研究，過深的巖漿洋會導致地幔中Ni、Co、Cr、V等元素的含量與觀測不符.雖然這個問題暫時還沒有答案，但最近的地月系統模擬似乎給出了新的思路：Lock和Stewart(2019)提出，如果大撞擊發生在高速旋轉的原始地球，則當時地球內部的壓力會比現在小一半左右，即撞擊后已經融合成現今大小的地球核幔邊界的壓力大約只有60～70 GPa左右，這正好與前人用現今地幔成分判據約束的巖漿洋底部壓力大致相符.

最后，值得注意的是，圖7中與五角星代表的“最優解”接近的巖漿洋應該在核幔邊界溫度為4560 K左右的時候，此溫度略高于巖漿洋的60%結晶溫度(Monteux et al., 2016).按照前人提出的地核熱演化模型(Liu et al., 2020; O′Rourke et al., 2017)，這似乎表明地核在42億年以前就已經通過核幔平衡獲得了現今的成分，此時地核中的鎂含量僅為0.3 wt%左右.因此，前人提出的通過硅析出或者鎂析出來解釋35億年以前的地球早期磁場就存在疑問了(Hirose et al., 2017; O′Rourke and Stevenson, 2016)，此時地核中生熱元素(U、Th、K等)的作用似乎又需要重新加以考慮(Chidester et al., 2022).

4 結語

核幔間的氧交換和地核中的氧含量是解譯地球演化歷史的關鍵線索，其底層物理化學機制在于核幔間元素的平衡.我們利用第一性原理模擬與機器學習技術，構建了超大尺度的兩相共存模擬，得到了核幔邊界壓力下硅酸鹽熔體和液態鐵平衡后的兩相組成，并結合實驗數據建立了新的平衡模型，這使得核幔間元素交換的研究可以更可靠地推至現今核幔邊界壓力.

運用本研究提出的平衡模型，我們深入探討了與核幔氧交換有關的熱點問題.我們發現只有在巖漿洋底部壓力達到80 GPa以上時，核幔平衡才能造就與觀測相符的現今地幔與地核主要元素的組成，因此大撞擊和超深巖漿洋事件對于解決核幔氧交換過程的爭議非常重要.

核幔元素交換問題非常復雜，本研究對約束此問題所做的努力顯然仍然不足.在模擬技術方面，我們雖然實現了超大體系的第一性原理分子動力學模擬，該方法仍有很多技術難題值得探索，比如地幔深處的含鐵硅酸鹽、氧化物礦物和熔體存在自旋轉變效應，這在我們的模擬中還沒有考慮，由于核幔元素交換實際上是各元素在不同相中化學勢的微小差別造成的，這些效應或許會顯著改變核幔間的平衡和元素分配.在平衡模型方面，我們目前的模型有多處做了簡化處理，比如平衡常數隨溫壓的近似變化關系、暫時未考慮硅酸鹽相中各組分的活度系數、陰陽離子價態變化、以及地核中其他更多元素的相互作用參數等.這些問題都將是我們今后研究的主要改進方向.

致謝感謝中國科學院地質與地球物理研究所超級計算中心、國家計算網格廣州中心以及北京超級云計算中心提供的支持與服務.感謝兩位審稿人和編輯對本文提出的寶貴意見.