章文博，張南，李昊遠

北京大學地球與空間科學學院，北京 100871

0 引言

月幔被認為由大碰撞后廣泛熔融的月球巖漿洋(Lunar Magma Ocean,LMO)結(jié)晶形成.巖漿洋逐漸冷卻凝固，結(jié)晶出不同的礦物，使得月幔自下而上形成了橄欖石、單斜輝石、斜方輝石的化學分層結(jié)構(gòu).最后階段，月幔會形成致密的富含鈦鐵礦等含鐵礦物的堆晶層 (以下簡稱鈦鐵礦層)，同時具有不相容性的KREEP也會逐漸在結(jié)晶的殘余液體中富集(Snyder et al.,1992;Elardo et al.,2011;Elkins-Tanton et al.,2011;Lin et al.,2017,2020;Rapp and Draper,2018).由于鈦鐵礦層密度遠大于下覆的月幔橄欖巖密度，因此堆晶層會發(fā)生重力失穩(wěn)下墜，此過程被稱為月幔倒轉(zhuǎn)(Lunar Cumulate Mantle Overturn)，許多天體演化的過程中都會經(jīng)歷類似的階段，比如地球和火星(Hess and Parmentier,1995;Parmentier et al.,2002).

巖漿洋結(jié)晶形成固態(tài)橄欖巖月幔和斜長巖月殼，之后月球繼續(xù)經(jīng)歷地質(zhì)演化.一方面，月表早期經(jīng)歷了頻繁的撞擊事件導致月表部分熔融，進而引發(fā)了許多巖漿事件(Norman,2009;Ding and Zhang,2018).另一方面，月幔部分熔融所產(chǎn)生的巖漿在表面噴發(fā)形成了月海玄武巖 (主要填充在地勢較低的區(qū)域，這些區(qū)域因為呈現(xiàn)暗色故被稱為月海)、月球玻璃等巖石(Delano,1986;Shearer et al.,2006).此外，盡管還有爭議(Mighani et al.,2020;Tarduno et al.,2021)，月球古強磁場的觀測建議月球內(nèi)部活動至少持續(xù)到1.5 Ga(Mighani et al.,2020)，在這之后月球內(nèi)部活動仍然持續(xù)，但強度已經(jīng)迅速衰減，其活動證據(jù)為月球深震所記錄(Khan et al.,2014;Weiss and Tikoo,2014).

月幔演化主要的地質(zhì)和動力學約束來自于月海玄武巖.月海玄武巖具有高FeO，低Al2O3，部分富含TiO2的巖石化學特征，根據(jù)TiO2的含量，被分為極低鈦(Very-low-Ti/VLT basalts,<1 wt.% TiO2)、低鈦(Low-Ti basalts,1～5 wt.% TiO2)以及高鈦玄武巖 (High-Ti basalts,>5 wt.% TiO2).空間上月海玄武巖的面積大約占月表面積的17%，主要分布在近地側(cè) (Head,1976).噴發(fā)時間主要集中于3.9～3.1 Ga(Nyquist et al.,2001;Hiesinger et al.,2011)，最年輕的月海玄武巖樣品年齡為1.96 Ga(圖1，Che et al.,2021).月海玄武巖的形成機制主要與月幔倒轉(zhuǎn)后富鈦物質(zhì)的上涌有關(guān)(Hess and Parmentier,1995;Zhong et al.,2000;Zhang et al.,2013a;Li et al.,2019).致密的鈦鐵礦層可以下沉進入月幔，甚至到達月球的核幔邊界(Core Mantle Boundary,CMB).隨著鈦鐵礦層中攜帶的產(chǎn)熱元素的衰變，發(fā)生倒轉(zhuǎn)的鈦鐵礦層會隨時間逐漸積累熱量并升溫.一方面，鈦鐵礦層加熱過程會影響月球核幔邊界的熱流，從而激發(fā)月核的能量變化，進而參與月球發(fā)電機的演化(Christensen et al.,2009;Weiss and Tikoo,2014)；另一方面，高溫的鈦鐵礦層會產(chǎn)生熱不穩(wěn)定性，在某些條件下會形成上涌的熱柱.在熱的鈦鐵礦物質(zhì)上升過程中，月幔也會被加熱熔融，鈦鐵礦物質(zhì)與月幔發(fā)生熔融與同化混染形成巖漿噴發(fā)到月表，可以解釋月海玄武巖和月球玻璃的形成機制(Elkins-Tanton et al.,2003;Grove and Krawczynski,2009;Barr and Grove,2013;Brown and Grove,2015).

圖1 Apollo玄武巖樣品,高地巖石樣品,玄武質(zhì)和斜長質(zhì)月球隕石樣品的年齡的總結(jié).定年方法在圖例中給出月海玄武巖的年齡主要來自ApolloBasalt DB_V2(Cone et al.,2020)以及Snape等(2019).玄武質(zhì)隕石包括Dhofar 287,NWA 032,NWA 773,KAL 009,MIL 13317,Y 793169,A 881757,MET 01210,MIL 05035,EET 96008,LAP 02205,NEA 003A和NWA 4734.年齡來自于Merle等(2020)和文中的參考文獻.月球玻璃包括Apollo 15綠色玻璃(Shih et al.,2001),Apollo 17橙色玻璃(Shih et al.,2001)和Apollo 17 VLT(Heiken et al.,1991).高地巖石年齡來自于Elkins-Tanton等(2011),Carlson等(2014),Borg等(2015)以及文中的參考文獻.高地隕石包括 Y 86032(Nyquist et al.,2006)和Dhofar 025(Leont′eva et al.,2005).CSFD (Crater Size-Frequency Distribution)方法的結(jié)果來自于Hiesinger等(2011).1.96 Ga為嫦娥五號著陸區(qū)玄武巖樣品U-Pb年齡(Che et al.,2021).Fig.1 Compilation of ages from Apollo mare basalts,highland rocks,basaltic and anorthositic meteorites.Analytical methods are shown Ages of mare basalts are mainly from ApolloBasalt DB_V2 (Cone et al.,2020)with the addition of Snape et al.(2019).Basaltic meteorites include Dhofar 287,NWA 032,NWA 773,KAL 009,MIL 13317,Y 793169,A 881757,MET 01210,MIL 05035,EET 96008,LAP 02205,NEA 003A and NWA 4734.Their ages are from Merle et al.(2020)and reference therein.Lunar glasses include Apollo 15 green glass (Shih et al.,2001),Apollo 17 orange glass (Shih et al.,2001)and Apollo 17 VLT (Heiken et al.,1991).Summary of highlands rock ages are from Elkins-Tanton et al.(2011),Carlson et al.(2014),Borg et al.(2015)and reference therein.Highland meteorites consist of Y 86032 (Nyquist et al.,2006)and Dhofar 025 (Leont′eva et al.,2005).The result of CSFD (Crater Size-Frequency Distribution)is from Hiesinger et al.(2011).1.96 Ga is the U-Pb age of mare basalts from Chang-E′5 landing site (Che et al.,2021).

除了月海玄武巖的分布以外，月球表面KREEP的空間分布也具有不對稱性.遙感觀測表明，月表K、Th、U等元素主要在風暴洋克里普巖地體(Procellarum KREEP Terrane,PKT)富集，而該區(qū)域與月海玄武巖的分布大致重合，因此兩者的空間分布很可能存在相關(guān)的作用機制(Lucey et al.,1998;Prettyman et al.,2006).

前人對月海玄武巖形成的內(nèi)動力過程研究主要集中于對鈦鐵礦層下沉的月幔倒轉(zhuǎn)以及下沉鈦鐵礦層受熱發(fā)生上涌兩個主要過程分別的數(shù)值模擬和參數(shù)探索，而將上述兩個過程統(tǒng)一起的數(shù)值模擬尚未被探索研究.

鈦鐵礦層的倒轉(zhuǎn)過程容易受到該層厚度和密度的影響(Parmentier et al.,2002;Li et al.,2019).這兩種性質(zhì)都依賴于巖漿洋的凝固過程.巖石學實驗和計算模擬實驗均給出了鈦鐵礦層的體積和結(jié)晶礦物組合的約束，前者對應(yīng)于鈦鐵礦層厚度，后者對應(yīng)于鈦鐵礦層密度.然而，由于這些實驗對巖漿洋的初始化學組成和深度有不同的假設(shè)，它們得出的鈦鐵礦層的厚度和密度會有一定程度的變化(Snyder et al.,1992;Elardo et al.,2011;Elkins-Tanton et al.,2011;Lin et al.,2017,2020;Rapp and Draper,2018).

內(nèi)部生熱在模型演化中也起著至關(guān)重要的作用.235U,238U,232Th和40K的放射性衰變是熱化學演化的主要熱源.同時，這些元素都具有不相容性，所以會在月球巖漿洋結(jié)晶的殘余液體中富集，也就是KREEP.盡管現(xiàn)在關(guān)于全硅酸鹽月球(Bulk Silicate Moon,BSM)的產(chǎn)熱元素豐度方面仍然存在分歧，但月球上U和Th元素的豐度被認為更可能接近于地球的豐度，而月球相對于地球更可能虧損鉀元素(Taylor,1982;Taylor and Wieczorek,2014).

前人的研究揭示了鈦鐵礦層的流變學性質(zhì)也對理解月球的熱化學演化有必不可少的幫助(Zhang et al.,2013a,b,2017;Dygert et al.,2016;Li et al.,2019).Zhang等(2013a)的數(shù)值模擬結(jié)果表明，如果鈦鐵礦層具有更高的活化能，那么被產(chǎn)熱元素加熱后，其黏度會更小，從而增強倒轉(zhuǎn)后的鈦鐵礦層內(nèi)部對流，反過來抑制熱量積累.因此，倒轉(zhuǎn)的鈦鐵礦層會更難產(chǎn)生熱不穩(wěn)定性，也就是說更難形成上涌熱柱.Zhang等(2017)考慮了鈦鐵礦的流變性質(zhì)，發(fā)現(xiàn)鈦鐵礦的軟化效應(yīng)會產(chǎn)生類似于高活化能的效果，這也不利于上涌熱柱的形成.Li等(2019)將化學組分和溫度相關(guān)的流變性質(zhì)應(yīng)用于月幔倒轉(zhuǎn)的數(shù)值模擬過程中，其結(jié)果表明鈦鐵礦對鈦鐵礦層黏度的減少會有利于月幔倒轉(zhuǎn)過程的進行.

本文的工作主要是同時結(jié)合了月幔倒轉(zhuǎn)和鈦鐵礦層上涌的過程，并探索鈦鐵礦層厚度、內(nèi)部生熱、瑞利數(shù)、成分和溫度相關(guān)的流變學的影響，進一步討論月海玄武巖的形成機制和相應(yīng)的生熱元素等的參數(shù)空間.

1 研究方法

1.1 控制方程與有限元數(shù)值模擬

月球熱化學演化的實現(xiàn)采用了三維球面幾何數(shù)值模型，模型假設(shè)無限大的普郎特數(shù)(Prandtl Number)，使用了布辛涅斯克近似(Boussinesq Approximation)，化學成分以及溫度相關(guān)的流變學和參數(shù)化的月核演化方程(Zhong et al.,2008;Zhang et al.,2013a,b,2017;Li et al.,2019).

無量綱形式的模型控制方程為：

(1)

(2)

(3)

(4)

(5)

(6)

其中α,ρ0,g,ΔT,R0,Rc,η0,κ和ρIBC分別是熱膨脹系數(shù)、月幔密度、月表重力加速度 (重力加速度隨深度增加而減小)、上下邊界的初始溫度差、月球半徑、月核半徑、參考黏度、熱擴散系數(shù)以及鈦鐵礦層密度.Ra和Rb可以用來定義浮力數(shù)B：

(7)

利用阿倫尼烏斯定律(Arrhenius Law)，我們計算鈦鐵礦層和月幔的混合密度：

(8)

其中ηIBC,E,R分別是鈦鐵礦在參考溫度下的黏度、有效活化能以及氣體常數(shù).

我們使用有限元計算軟件CitcomS來數(shù)值求解三維球殼中的不可壓縮流方程(1)—(4)，用來模擬月幔固結(jié)后內(nèi)部物質(zhì)的重力倒轉(zhuǎn)和加熱產(chǎn)生熔體上涌噴發(fā)的過程.球殼從下到上分別為核幔邊界、月幔層、鈦鐵礦層、月殼層以及月表，三維球殼被均分成12個區(qū)域，每個區(qū)域有65×65×65結(jié)點.為了提高分辨率，鈦鐵礦層的結(jié)點在徑向上被進一步細分，根據(jù)其厚度，結(jié)點數(shù)量從5～15不等.使用自由滑動條件作為上下邊界的力學邊界條件.依照Zhang等(2013a)的處理方式，月表溫度的無量綱值恒為1，CMB的初始無量綱溫度為1，之后隨著月核與月幔之間的能量守恒關(guān)系演化.

1.2 參數(shù)空間

在我們的模型中，有四個關(guān)鍵參數(shù)需要被探索，分別是：鈦鐵礦層厚度、內(nèi)部生熱、參考黏度 (瑞利數(shù))以及鈦鐵礦層黏度.其他參數(shù)在表1中給出.

表1 模型參數(shù)Table 1 Model parameters

我們的模型研究了50～150 km的鈦鐵礦層厚度.對于1000 km深度的LMO深度，通過質(zhì)量守恒和完全分離結(jié)晶的相平衡關(guān)系得到的鈦鐵礦層厚度范圍是18～50 km，對應(yīng)的鈦鐵礦層中鈦鐵礦體積占比為16 vol.%到3.6 vol.%(Snyder et al.,1992;Elardo et al.,2011;Elkins-Tanton et al.,2011;Lin et al.,2017,2020;Rapp and Draper,2018;Li et al.,2019).然而，考慮到在LMO完全固結(jié)之前月幔倒轉(zhuǎn)就已經(jīng)啟動了，因此我們引入固定鈦鐵礦總量但是對體積進行一定程度稀釋的鈦鐵礦層.基于上述文獻提到的鈦鐵礦層厚度以及鈦鐵礦含量的上下界，我們考慮厚度和鈦鐵礦含量的算數(shù)平均值，即20 km的鈦鐵礦層厚度和11 vol.%的鈦鐵礦含量，并對其進行稀釋 (圖2a).

圖2 (a)鈦鐵礦含量與鈦鐵礦層厚度的關(guān)系,修改自 Li等(2019).淺藍色實線代表50 km 和3.6 vol.%的下界.淺藍色虛線代表18 km和16 vol.%的上界.深藍色實線為算數(shù)平均值20 km和11 vol.%.(b)月球內(nèi)部生熱的四套方案.兩條實線分別為Turcotte和Schubert(2002)、McDonough和Sun(1995)的結(jié)果,兩條虛線為對應(yīng)的忽略鉀元素的結(jié)果.Taylor′s Moon的產(chǎn)熱量用黑色點劃線表示，與Turcotte和Schubert (2002)無鉀方案相當，即紅色虛線Fig.2 (a)The derived relation between the ilmenite volume fraction and IBC thickness,from Li et al.(2019).Solid light blue line is the lower bound with 50 km and 3.6 vol.%.Dashed light blue line is the upper bound with 18 km and 16 vol.%.Solid dark blue line is the average value with 20 km and 11 vol.%.(b)Four schemes of lunar internal heating.Solid lines are original estimates from Turcotte and Schubert (2002)and McDonough and Sun (1995),respectively.Dash lines are corresponding K-free results.Dash-dotted black line is heat-producing rate of Taylor′s Moon,which is equivalent to the K-free scheme of Turcotte and Schubert (2002),that is,the dashed red line

正如之前提到的，鈾、釷元素作為難熔元素，其在BSM的含量與全硅酸鹽地球(Bulk Silicate Earth,BSE)的含量接近，但是作為揮發(fā)性元素的鉀元素在月球中的含量則有很大的不確定度(Taylor and Wieczorek,2014).因此，對鈾、釷元素，我們采用地球自然豐度.但是目前對全硅酸鹽地球的放射性同位素含量估計也存在一些分歧，于是我們考慮兩套方案：鈾、釷元素豐度的下界來自于McDonough和Sun(1995)的地幔巖模型，而上界則來自于Turcotte和Schubert(2002)的地球物理模型(表2).相對于地球，月球很可能高度虧損鉀元素，因此我們?nèi)サ羯鲜鰞蓚€模型的鉀元素豐度，從而導出兩個無鉀的方案.為了便于討論內(nèi)部生熱的效應(yīng)，我們將以上四套產(chǎn)熱方案統(tǒng)一到鈾、釷元素豐度分別為20.3 ng·g-1和79.5 ng·g-1的情況，相應(yīng)的K/U比值分別為0，8500，11800和23500(圖2b).前人的數(shù)值和化學模型所使用的豐度含量大多落在K/U值為0～8500的范圍內(nèi)(Zhong et al.,2000;Zhang et al.,2013a,b,2017;Laneuville et al.,2018;Li et al.,2019).此外，我們也考慮了月球相對地球更富集鈾、釷元素的模型，即 Taylor′s Moon(Taylor et al.,2006)，發(fā)現(xiàn)其產(chǎn)熱量與K/U=8500的產(chǎn)熱方案相當(圖2b).被鈦鐵礦攜帶的KREEP含量目前仍然缺少實驗約束，因此我們采用了Zhong等(2000)和Zhang等(2013a)的假設(shè)，即模型開始時所有的KREEP與鈦鐵礦層完全混合.產(chǎn)熱元素的衰變常數(shù)依據(jù)Turcotte和Schubert(2002).

表2 U,Th和K的預(yù)計濃度Table 2 Estimated concentrations of U,Th and K

我們將參考溫度Tref設(shè)為1300 ℃，對應(yīng)的橄欖巖黏度η0范圍為5×1018到1×1021Pa·s(Hirth and Kohlstedt,2004).考慮1×1020Pa·s,5×1020Pa·s和1×1021Pa·s的參考黏度，對應(yīng)的熱瑞利數(shù)為7.6×106,1.5×106和7.6×105.鈦鐵礦層的黏度由鈦鐵礦含量、含水量以及物理狀態(tài)比如熔融和溫度決定 (Dygert et al.,2016).取鈦鐵礦層相對于橄欖巖的黏度為1×10-1到1×10-2.此外，考慮橄欖巖和鈦鐵礦層具有相同的活化能100 kJ·mol-1，雖然鈦鐵礦和橄欖石的流變學實驗得到的活化能更大(Karato and Wu,1993;Hirth and Kohlstedt,2004;Dygert et al.,2016)，但那是在非線性流變的情況下得到的值，Christensen(1984)認為在假設(shè)應(yīng)變速率不變的情況下，將實驗得到的活化能乘以0.3，可以將非線性流變的效果轉(zhuǎn)化為線性流變的效果，約為100 kJ·mol-1，并由此定義出有效黏度.

1.3 初始條件

模型開始的初始年齡為 4.42 Ga，盡管最近的研究表明LMO完全結(jié)晶的年齡可能更晚，但是我們的模型仍然將月球斜長巖樣品的平均年齡作為LMO結(jié)晶末期年齡，即4.42 Ga(Elkins-Tanton et al.,2011;Borg et al.,2015,2019;Maurice et al.,2020).我們的模型采用了和Zhang等(2013a)一樣的橄欖巖固相線溫度作為初始溫度剖面(圖3c)，模型初始密度分布在圖3d中給出.算例的命名方式會給出鈦鐵礦層的厚度、相對黏度、參考黏度以及K/U比值.比如對于算例“D50v1e-1V1e20K0”，即代表初始的鈦鐵礦層厚度為50 km，相對黏度為1×10-1，參考黏度為1×1020Pa·s，K/U比值為0.

圖3 (a)鈦鐵礦層厚度為50 km和75 km的算例初始黏度剖面.(b)鈦鐵礦層厚度為100 km和150 km的算例初始黏度剖面.(c)初始溫度剖面以及四個算例即將產(chǎn)生上升熱柱時的溫度剖面.(d)四種鈦鐵礦層厚度的初始密度剖面Fig.3 (a)Initial viscosity profiles for cases of 50 km and 75 km.(b)Initial viscosity profiles for cases of 100 km and 150 km.(c)Initial temperature profile and temperature profiles in the beginning of upwelling for four cases.(d)Initial density profiles for four sets of IBC thickness

1.4 玄武巖上涌時間定義

為了研究發(fā)生月幔倒轉(zhuǎn)后的內(nèi)部演化，我們計算了演化過程中，在深度為450 km處的球面上，TiO2含量大于0.46 wt.%的結(jié)點占比，我們認為結(jié)點占比超過15%是成功的上涌熱柱的標志，基于此，成功上涌的時間被定義為：

tupwelling=t[ratio=15%]-4.42 Ga,

(9)

我們使用該簡化準則來判斷是否形成了上涌熱柱，理由如下：(1)月海玄武巖來源于更為分異的巖漿，而與玄武巖共存的月球玻璃則被認為代表了更為原始的巖漿成分(Shearer et al.,2006；Xiao et al.,2021).本模型不考慮復(fù)雜的化學演化，因此將鎂鐵質(zhì)月球玻璃的組分作為玄武巖的替代指標.月球玻璃的多相飽和壓強(Multiple Saturation Pressures)可以通過高溫高壓實驗確定，前人的實驗表明，月球玻璃的最大深度約為450 km，其對應(yīng)的熔體的TiO2含量為0.46 wt.%(圖4，Delano,1986;Elkins-Tanton et al.,2003;Grove and Krawczynski,2009).因此，我們采用0.46 wt.%作為模型計算的玄武巖TiO2含量.(2)對月表玄武巖的觀測表明，月海玄武巖的面積占月表面積的16%～17%，并且其噴發(fā)年齡主要分布于3.8～3.3 Ga之間(Head,1976;Whitten and Head,2015).為了模擬月海玄武巖母巖漿的生成，我們將月海玄武巖所占面積的比例大致對應(yīng)于在源區(qū)450 km深度時的面積比例，即在450 km深度達到(1)中的TiO2含量的結(jié)點比例.

圖4 綠色玻璃的密度、多相飽和壓強以及相應(yīng)TiO2含量，圖修改自 Grove和Krawczynski (2009)，TiO2含量來自Delano(1986)Fig.4 Green glasses′ density,multiple saturation pressure from Grove and Krawczynski (2009)and contents of TiO2 from Delano (1986)

2 結(jié)果

2.1 鈦鐵礦層厚度的影響

首先，我們研究了鈦鐵礦層厚度對上涌時間的影響.當參考黏度為1×1020Pa·s,K/U比為0,鈦鐵礦層相對黏度為1×10-1時，分別指定鈦鐵礦層厚度為50 km、75 km、100 km和150 km(圖5a).

圖5 不同鈦鐵礦層厚度，K/U 比值，參考黏度下的上涌時間圓圈、叉號以及三角形符號分別代表成功上涌的算例、無法形成上涌熱柱的算例以及根據(jù)規(guī)律推測的無法形成上涌的參數(shù)空間.四條豎線分別代表圖2b中的四種生熱方案.(a),(b),(c)分別為不同參考黏度下的結(jié)果.(b)和(c)中的粗虛線是能否成功形成上升熱柱的分界線.Fig.5 Upwelling time with respect to the IBC thickness,K/U ratio,reference viscosity Circles,crosses and triangles are cases of successful upwelling,failed upwelling and deduced failed upwelling according to upwelling patterns,respectively.Four vertical lines correspond to four schemes in Fig.2b.(a),(b),(c)are results of different reference viscosity.Thick dashed lines in (b)and (c)are dividing lines between successful upwelling and that of failed upwelling.

標準算例為D50v1e-1V1e20K0.熱化學演化包括三個階段.開始階段，致密的冷的鈦鐵礦層下沉(圖6a,6b).鈦鐵礦層將一直下沉直到核幔邊界深度，最終經(jīng)過約200 Ma后，48.5 vol.%的鈦鐵礦層會進入下月幔.第二個階段，在底部聚集的鈦鐵礦層會停留約100 Ma，期間鈦鐵礦層所攜帶的產(chǎn)熱元素會不斷衰變并加熱鈦鐵礦層(圖6c，6d).在約4.12 Ga時，被加熱的鈦鐵礦層產(chǎn)生了足夠的浮力，開始形成上涌熱柱(圖6e，6f).在接下來的100 Ma間，底部的鈦鐵礦層會不斷地被熱能驅(qū)動，形成很強的上涌熱柱，熱柱的三維形態(tài)經(jīng)過三維球諧函數(shù)分解，主要的球諧階數(shù)為3階(圖6g，6h).

圖6 標準算例D50v1e-1V1e20K0的三維熱化學結(jié)構(gòu)以及相應(yīng)的二維溫度剖面(a)和(b)是倒轉(zhuǎn)階段.(c)和(d)是加熱階段.(e)和(f)是上涌的啟動階段.(g)和(h)是上涌階段.三維結(jié)構(gòu)中，透明的綠色等值面代表化學組分為0.14的面.紅色和藍色等值面分別代表溫度比球面平均高250 ℃，低250 ℃的面.白色的球面代表月核.在(h)中,箭頭代表流場速度.三維和二維圖中的藍色圓代表450 km深度.Fig.6 3-D thermochemical structures and corresponding 2-D cross sections for the standard case D50v1e-1V1e20K0(a)and (b)are overturned phase.(c)and (d)are heating phase.(e)and (f)are onset of upwelling.(g)and (h)are upwelling phase.In the 3-D view,transparent green iso-surfaces represent the chemical composition at a contour level of 0.14.Red and blue iso-surfaces represent the residue temperature with the contour levels at 250 ℃,-250 ℃.White sphere stands for the lunar core.In (h),the arrows show the flow velocity.Blue circles in both 3-D and 2-D graphs stand for the 450-km depth.

與D50v1e-1V1e20K0相比，更厚的鈦鐵礦層，即算例D75v1e-1V1e20K0,D100v1e-1V1e20K0,D150v1e-1V1e20K0有著相似的演化模式.但是，也存在一些值得關(guān)注的區(qū)別.由上涌時間的變化可以看出，其明顯與鈦鐵礦層厚度負相關(guān).隨著鈦鐵礦層厚度從50增大到150 km，上涌時間從387 Ma單調(diào)縮短到354 Ma.對上涌時間的影響來源于月幔倒轉(zhuǎn)、加熱和上涌過程各自的貢獻.一方面，更厚的鈦鐵礦層有利于月幔倒轉(zhuǎn)過程的進行，這不僅會增加進入下月幔鈦鐵礦層的比例，形成更厚的底部鈦鐵礦層，還會縮短發(fā)生月幔倒轉(zhuǎn)的時間尺度(圖7a).另一方面，更厚的鈦鐵礦層會導致更小的化學密度，而這有利于底部鈦鐵礦層更快地產(chǎn)生上涌熱柱.

2.2 內(nèi)部生熱的影響

其次，我們探索了內(nèi)部生熱的影響.與2.1節(jié)類似，固定鈦鐵礦層厚度為50 km，鈦鐵礦層相對黏度為1×10-1，參考黏度為1×1020Pa·s，分別指定K/U比值為0、8500、11800和23500(圖5a).

K/U比值的增加等效于內(nèi)部生熱的增加，因此有利于形成上涌熱柱.上涌時間從387 Ma縮短到184 Ma.除了演化時間以外，其他指標如發(fā)生倒轉(zhuǎn)的鈦鐵礦層比例則幾乎不受內(nèi)部生熱影響.因此，增加內(nèi)部生熱僅僅會加速加熱階段，而對其他階段的影響很小(圖7b，圖8).

圖7 下月幔中鈦鐵礦層比例隨時間的演化垂直虛線代表上涌時間點，下同.(a)代表不同鈦鐵礦層厚度的四個算例.(b)代表不同內(nèi)部生熱的四個算例.Fig.7 Evolutions of overturned IBC fraction in the lower lunar mantleVertical dashed lines stand for the upwelling time and are the same as following figures.(a)is for four cases of different IBC thickness.(b)is for four cases of different internal heating.

圖8 內(nèi)部生熱對平均溫度的影響的演化圖Fig.8 Evolutions of the azimuthally averaged temperature with respect to internal heating

2.3 月幔參考黏度的影響

參考黏度/瑞利數(shù)對熱化學演化也有顯著的影響.考慮兩組算例，第一組算例為D150v1e-1V1e20K0,D150v1e-1V5e20K0和D150v1e-1V1e21K0；第二組算例為D50v1e-1V1e20K23500,D50v1e-1V5e20K23500和D50v1e-1V1e21K23500.

對第一組算例，可以看出上涌時間與參考黏度呈正相關(guān)關(guān)系.從1×1020Pa·s到1×1021Pa·s，參考黏度增加了一個數(shù)量級，相應(yīng)的上涌時間增加了一倍，即從354 Ma增加到了793 Ma(表3,圖9a、b、c).因此，增加參考黏度會極大推遲上涌熱柱的形成.

更重要的是，第二組算例證明了增加參考黏度甚至會阻止熱柱形成.由于K23500是我們考慮的最大產(chǎn)熱方案，但僅有算例D50v1e-1V1e20K23500在184 Ma后形成了上升熱柱(表3,圖9d、e、f).實際上，增加參考黏度等效于減小瑞利數(shù)，而瑞利數(shù)是對流系統(tǒng)活躍程度的表征.除了上涌時間以外，倒轉(zhuǎn)的鈦鐵礦層比例也在很大程度上受到參考黏度的影響.對于以上考慮的兩組算例，黏度從1×1020Pa·s,5×1020Pa·s到1×1021Pa·s變化時，倒轉(zhuǎn)的鈦鐵礦層比例也減小為原來的50%(表3,圖9).最大的倒轉(zhuǎn)比例發(fā)生在1×1020Pa·s的算例中，分別是69.3%和48.7%，而最小的比例則發(fā)生在1×1021Pa·s的算例中，對應(yīng)比例分別減小到35.5%和22.0%(表3，圖9).因此，當瑞利數(shù)減小時上涌熱柱更快形成的結(jié)果是符合預(yù)期的.

圖9 參考黏度對平均密度的影響的演化圖Fig.9 Evolutions of azimuthally averaged density with respect to reference viscosity

2.4 鈦鐵礦層黏度的影響

最后，我們分析了鈦鐵礦層黏度的影響效果，考慮了四個算例：D75v1e-1V5e20K11800,D75v1e-2V5e20K11800,D100v1e-1V5e20K8500和D100v1e-2V5e20K8500.

與鈦鐵礦層厚度、內(nèi)部生熱以及參考黏度不同，鈦鐵礦層黏度對于上涌時間的影響更為復(fù)雜.對于D75的兩個算例，鈦鐵礦層黏度從1×10-1Pa·s到1×10-2Pa·s的減小會導致上涌失敗(表3).相反，在D100的兩個算例，當鈦鐵礦層黏度同樣減小一個數(shù)量級時，上涌時間從449 Ma減小到325 Ma.

表3 關(guān)鍵輸入輸出參數(shù)Table 3 Key input and output parameters

解釋上述效應(yīng)需要考慮兩個因素，分別是底部鈦鐵礦層的厚度以及化學組分的混合效應(yīng).首先，更小的鈦鐵礦層黏度有利于倒轉(zhuǎn)過程的進行，因為鈦鐵礦層黏度越小，在受重力影響進入下月幔的過程中會受到更小的阻力，使得倒轉(zhuǎn)過程更早啟動，并且有更多的鈦鐵礦層會發(fā)生倒轉(zhuǎn)(圖10).對于更厚的底部鈦鐵礦層，熱量會更容易聚集起來，從而有利于上涌熱柱的形成(圖3c).其次，化學組分的混合效應(yīng)也對上涌機制有著重要的影響.但鈦鐵礦層的相對黏度為1×10-1時，我們注意到鈦鐵礦層和月幔發(fā)生了更高程度的混合，而相對黏度為1×10-2時，鈦鐵礦層與月幔混合程度則更低(圖11b，d，f，h).發(fā)生混合的鈦鐵礦層比純的鈦鐵礦層密度更小，因此在加熱過程中更容易獲得足夠的浮力.所以，在某種程度上，鈦鐵礦層混合越充分，越容易形成上涌熱柱.

圖10 鈦鐵礦層黏度對進入下月幔的鈦鐵礦層比例的影響的演化圖Fig.10 Evolutions of overturned IBC fraction in the lower lunar mantle with respect to IBC viscosity

圖11 鈦鐵礦層黏度的影響中的4個算例在上涌之處的三維形態(tài)結(jié)構(gòu)圖和對應(yīng)二維化學組分圖三維圖(a),(c)和(f)的綠色透明等值面代表化學組分為0.21的面，在(e)中為0.27的等值面.其他均與圖6相同.Fig.11 3-D thermochemical structures and corresponding 2-D chemical cross sections for 4 cases in the beginning of upwelling to test the effect of IBC viscosityIn the 3-D view,transparent green iso-surfaces in (a),(c)and (f)represent the chemical composition at a contour level of 0.21 and that in (e)is at a contour level of 0.27.Other color schemes are the same as those in Fig.6.

總而言之，鈦鐵礦層黏度更小，會導致底部鈦鐵礦層厚度更大，但是化學混合程度更弱.這兩個效應(yīng)對于上涌時間是相互對抗的.對于D100的兩個算例，鈦鐵礦層黏度對于底部厚度的影響效果超過了化學混合的效應(yīng)(圖11e,f,g和h)，而對于D75的兩個算例，兩者效果的相對權(quán)重則發(fā)生了顛倒(圖11a,b,c和d).

3 討論

月海玄武巖主要集中在近地側(cè)區(qū)域，這種不對稱性需要長波的上涌，即主導球諧階數(shù)最小(Zhong et al.,2000;Zhang et al.,2013a,2017).基于當前已知的鈦鐵礦的黏度，很遺憾我們暫時無法產(chǎn)生一個包含倒轉(zhuǎn)和上涌過程的球諧一階結(jié)構(gòu)，最接近長波上涌的算例為D150v1e-1V1e21K0，其上升熱柱的形態(tài)經(jīng)過球諧函數(shù)分解主要為二階(表3,圖12).

圖12 最長波長的算例的三維形態(tài)結(jié)構(gòu)圖以及二維溫度剖面三維圖中綠色透明等值面代表化學組分為0.4的面.其余均與上同.Fig.12 The 3-D thermochemical structures and 2-D cross sections for the case of the largest upwelling wavelengthIn the 3-D view,transparent green iso-surfaces represent the chemical composition at a contour level of 0.4.Other schemes are the same as those above.

重力失穩(wěn)的解析解也可以通過瑞利-泰勒不穩(wěn)定分析得到(圖13a)，其結(jié)果表明一階的上涌要么需要更小的底部鈦鐵礦層厚度(至多200 km)和鈦鐵礦層黏度(至多10-3的相對黏度)，要么需要黏度更大的鈦鐵礦層(至少101的相對黏度).因此，D150v1e-1V1e21K0落在了一階上涌的參數(shù)空間外(圖13a,b).如果要將鈦鐵礦層黏度減小到10-3，那么根據(jù)之前的結(jié)果，很難形成上涌熱柱.而如果要提高鈦鐵礦層黏度，我們需要將其他因素納入流變學性質(zhì)的考慮中，最新的巖石力學實驗發(fā)現(xiàn)鈦鐵礦有較大的活化體積(Tokle et al.,2021)，這意味著隨著深度(壓強)增加，鈦鐵礦的黏度會大大增加，希望這能夠幫助我們產(chǎn)生更加長波的月幔熱柱上涌.

圖13 (a)通過瑞利-泰勒不穩(wěn)定分析得到的在不同鈦鐵礦層黏度和底部鈦鐵礦層厚度下發(fā)生上涌時的主導球諧階數(shù).(b)兩個最長波長算例在即將形成上升流時的黏度剖面Fig.13 (a)The predicted dominant spherical harmonic degree of overturn as a function of the IBC viscosity contrast and the thickness of bottom IBC.(b)Viscosity profiles in the beginning of upwelling for two cases of the largest upwelling wavelength

熔融對模型演化的影響主要體現(xiàn)在兩個方面，在力學效應(yīng)方面，鈦鐵礦層受熱發(fā)生部分熔融，產(chǎn)生熔體增加巖石格架的孔隙率，并會極大的降低鈦鐵礦層的黏度；在熱效應(yīng)方面，鈦鐵礦層熔融時，熔體以潛熱形式吸收熱量，從而影響模型的熱演化過程.巖石學實驗表明對于部分熔融的橄欖巖，其應(yīng)變速率與eα φ成正比，其中φ是部分熔融程度，α是部分熔融系數(shù)，根據(jù)形變機制不同，取值在21～45之間(Hirth and Kohlstedt,2004;Mei et al.,2002;Scott and Kohlstedt,2006)，據(jù)此約10%的部分熔融可以導致黏度減小一個數(shù)量級.但是以上估計是基于熔體不發(fā)生運移的情況，實際情況下，固體熔融產(chǎn)生熔體后，固體格架會發(fā)生改變，同時熔體與固體存在密度差，因此產(chǎn)生的液相熔體會立刻發(fā)生運移和分隔，這使得，想要定量研究熔融對黏度的影響，需要利用考慮液相熔體的雙相流模型(例如Zhang et al.,2020).前人主要考慮熔融對熱演化的影響(Laneuville et al.,2018;Li et al.,2019)，我們也嘗試了一個只考慮熔融對熱演化影響的算例，為D150v1e-1V1e20K0M2(表3)，參照Li等(2019)，只改變控制方程中的能量方程，其余不變：

(10)

(11)

圖14 熔融(D150v1e-1V1e20K0M2)對進入下月幔的鈦鐵礦層比例的影響的演化圖Fig.14 Evolutions of overturned IBC fraction in the lower lunar mantle with respect to melting (D150v1e-1V1e20K0M2)

目前現(xiàn)存的月球巖石樣品、隕石樣品以及隕石坑數(shù)量定年方法均認為月球火山作用是一個長期過程，歷時超過10億年(Nyquist et al.,2001;Hiesinger et al.,2011;Joy and Arai,2013).阿波羅任務(wù)獲得的月球玄武巖樣品年齡主要分布于3.1～3.9 Ga，而玄武質(zhì)的月球隕石樣品年齡則擴展至4.37 Ga(KAL 009)到2.7 Ga(NWA 032,773)(圖1)(Nyquist et al.,2001;Fagan et al.,2002;Terada et al.,2007;Joy and Arai,2013;Snape et al.,2019)，中國嫦娥5號任務(wù)的著陸區(qū)的玄武巖樣品年齡為1.96 Ga(Che et al.,2021).然而，通過隕石坑數(shù)量定年得到的玄武巖年齡會更為年輕，可達1.2 Ga(Hiesinger et al.,2011,Qian et al.,2021).月球火山作用的時間分布特征需要一個能夠長期驅(qū)動月幔熔融的有效機制.我們的模型可以得到上涌熱柱的幕式行為(圖15).從3.9 Ga到3.3 Ga，450 km深度到初始鈦鐵礦層底部之間的鈦鐵礦層比例三次達到了極值，表示在這期間有三次主要的上涌，從而可能形成三次月幔熔融和玄武巖噴發(fā)事件(Qin and Zhong,2014).但我們的模型產(chǎn)生的最后一次上涌為3.3 Ga，這個時間比觀測到最新的月海玄武巖(1.96 Ga)還有很大的時間差別，所以可能需要更多的熱源來驅(qū)動月球內(nèi)部的熱演化.此前的熱演化研究中，月核與月幔差異運動導致的耗散功率常常被忽略(Williams et al.,2001;Dwyer et al.,2011)，Dwyer等(2011)估計了月球早期核幔差異運動的耗散功率，得到在4.0 Ga左右其耗散功率的量級為1×1012W.但是他們在討論中認為可能高估了耗散功率，因此又定性地討論了核幔耦合作用下，月球演化早期差異運動摩擦產(chǎn)生的功率為3×1011W.通過對本文中的四套生熱方案的單位生熱率(圖2b)可以估計得到四套方案的總生熱功率分別為18×1011W，14×1011W，12×1011W，7×1011W.所以核幔摩擦這一熱源對月球內(nèi)部更長期的熱演化有一定作用，我們認為他們的估計還不夠精確，我們將在以后的工作中精確評估核幔差異運動的效果，來對比和生熱元素貢獻的大小.

圖15 2000 Ma間在450 km深和初始鈦鐵礦層層底部之間的鈦鐵礦層比例隨時間的演化圖橙色陰影區(qū)標記了主要的上涌階段.Fig.15 Evolutions of upwelling IBC fraction between the 450-km depth and the bottom of initial IBC layer for 2000 MaOrange shadow regions denotes major upwelling periods with peaks of IBC fraction.

通過改變生熱元素豐度和參考黏度，模型給出的上涌時間可以用來限定這兩者的參數(shù)空間.我們的模型用到的參考黏度范圍為1×1020Pa·s到1×1021Pa·s，對地球上地幔的實驗結(jié)果給出黏度的范圍為1×1018Pa·s到1×1021Pa·s(Hirth and Kohlstedt,2004)，而前人模型中考慮的參考黏度范圍為5×1019Pa·s到1×1022Pa·s(Zhang et al.,2013a;Laneuville et al.,2018;Qin et al.,2018;Li et al.,2019;Zhao et al.,2019).我們的結(jié)果表明當參考黏度為1×1021Pa·s時，只有位于粗虛線上方的參數(shù)空間才能發(fā)生上涌，進一步來說，如果要月球主要玄武巖來自鈦鐵礦層倒轉(zhuǎn)又加熱上升的機制，而且集中在3.9 Ga大規(guī)模噴發(fā)，需要極厚的初始鈦鐵礦厚度(>100 km)和較大的生熱元素豐度，至少K/U為8500，逼近地球生熱元素豐度(圖5c).而如果參考黏度達到1×1022Pa·s(Qin et al.,2018)，對我們模型改變參考黏度得到上涌時間進行外推，可知該參考黏度下無法在3.9 Ga產(chǎn)生大規(guī)模玄武巖巖漿，Qin等(2018)也認為他們基于月球赤道突起而約束得到的月幔黏度1×1022Pa·s實際上對黏度范圍1×1021～1×1022Pa·s沒有分辨率，因此我們認為月幔參考黏度最大應(yīng)為1×1021Pa·s.總而言之，我們的模型表明為了鈦鐵礦和KREEP能重新噴發(fā)，月幔的參考黏度上限應(yīng)在1×1021Pa·s這個量級，此時月球生熱元素的豐度至少為Taylor′s Moon的程度，接近McDonough和Sun(1995)的地球豐度；當月幔參考黏度更小時，允許月球的生熱元素豐度小于McDonough和Sun(1995)的地球豐度.

4 結(jié)論

總而言之，我們使用了三維的有限元模型來研究月幔倒轉(zhuǎn)和倒轉(zhuǎn)之后的熱化學演化得到了如下結(jié)論.我們發(fā)現(xiàn)內(nèi)部生熱對玄武巖產(chǎn)生的影響較大，對于相同的鈦鐵礦層厚度、黏度以及月幔參考黏度，內(nèi)部生熱越強，熱柱上涌產(chǎn)生玄武巖巖漿所用的時間越短，反之若內(nèi)部生熱越弱，則所用時間越長.

我們通過月球在3.9 Ga集中噴發(fā)月海玄武巖的觀測結(jié)論約束得到了月幔參考黏度的上限為1×1021Pa·s，此時對于較厚的初始鈦鐵礦層，月球生熱元素的豐度至少應(yīng)該達到Taylor′s Moon的程度，如果月球初始鈦鐵礦層較薄 (<100 km)，則需要更多的生熱元素豐度，即McDonough和Sun(1995)的地球豐度或更高.

我們考慮了熔融對月球內(nèi)部熱演化的影響，結(jié)果表明，如果鈦鐵礦層發(fā)生熔融，那么將要求月球內(nèi)部有更高的生熱元素豐度或者其他貢獻較大的熱源，才能滿足3.9 Ga大規(guī)模玄武巖噴發(fā)的觀測約束.

本文還得到了月海玄武巖周期性的噴發(fā)模式以及玄武巖的不對稱分布，但是與實際月海玄武巖的時間和空間分布還有一定差距，未來的工作將繼續(xù)探索參數(shù)空間來接近實際情況.

致謝本研究的資金來源于CNSA D020205.感謝2位匿名評審人的詳細意見.其中建議的文獻提高了本文的視角.本文的大型并行計算在北京大學高性能計算平臺上完成.本文通訊作者感謝本科時期王仁老師的言傳身教，他在科研工作上的耐心和仔細給學生時期的作者留下深刻的印象.