三維黏彈性模型在月球冷卻產生的熱應力計算中的應用*

金一民,陶莎,石耀霖

(中國科學院大學地球與行星科學學院 中國科學院計算地球動力學重點實驗室, 北京 100049) (2022年1月28日收稿; 2022年6月8日收修改稿)

月球沒有活躍的板塊運動,其冷卻產生的熱應力得以長時間積累,因此熱應力成為影響月球內部力學環境的重要因素之一。由于月球內部的冷卻速度大于表面,故熱應力的總體分布為徑向拉張和切向擠壓[1]。研究者對美國宇航月球勘測衛星(lunar reconnaissance orbiter)在月球表面觀測到的逆沖斷層進行了分析,推測它們是月球冷卻產生的切向擠壓應力的產物[2-3]。同時,由于阿波羅月震實驗(Apollo lunar seismic experiments package)記錄到的28次淺源月震中有8次落在被解釋為逆沖斷層的地形陡坎附近30 km以內[4],因此月球淺部的熱應力釋放可能是觸發淺源月震的原因之一。

除淺源月震,發生在半徑600～1 000 km區域內的深源月震也可能與月球冷卻產生的熱應力存在關聯。由于深源月震數據與月球潮汐表現出相同的周期性(27 d的月地軌道周期和206 d的太陽擾動周期)[5-6],因此普遍認為深源月震與月球的潮汐應力相關[7-8]。然而,張貝等[9]的計算表明,若假設月球為分層彈性介質,則月幔深部的潮汐應力遠未達到足以觸發月震的數量級。為解釋潮汐作用觸發深源月震的機制,張貝等[9]提出在深源月震多發區存在高壓孔隙流體的假設,認為孔隙流體分擔了大部分靜巖壓力,使得發生破裂所需的差應力減小。月球內部普遍存在的徑向拉張熱應力為月幔深部孔隙的發育提供了一種可能的解釋,將在第4節做詳細討論。

盡管月球冷卻產生的熱應力可能是月震活動的影響因素之一,但由于對月球各圈層的力學性質缺乏足夠了解,目前對月球熱應力積累的定量計算很少。我們使用一維參量化模型[10]計算了月球的冷卻和熱應力積累過程,計算結果表明月球冷卻產生的熱應力在淺部以切向擠壓為主,足以產生逆沖斷層并觸發淺源月震,而深部的差應力較小,因此熱應力不是觸發深源月震的直接原因。受計算方法的限制,我們的一維參量化模型采用彈性本構關系計算熱應力,而對大時間尺度下不可忽略的黏性效應只做了簡化定量計算。本研究使用三維黏彈性有限元模型計算月球熱應力的積累過程,通過對比實驗分析黏性參數對熱應力的影響,并進一步討論熱應力和月震之間的關系。

1 數值模型

1.1 控制方程

若將月球的冷卻過程用一階差分格式離散成有限個時間步,則每個時間步內的力學狀態可以近似看作平衡狀態。第n個時間步的平衡方程可表示為

(1)

(2)

本實驗采用Maxwell黏彈性模型模擬月殼和月幔的力學性質[11-12],其本構方程為

其中：Δεn為應變增量,Δtn為第n個時間步的大小,α為線膨脹系數,ΔTn為溫度增量,I為單位矩陣。C和Q分別表示彈性應力-應變張量和黏性應力-應變率張量。對于三維各向同性介質,有

Cijkl=μ(δikδjl+δilδjk)+λδijδkl,

(4a)

其中：λ和μ為Lamé常數,η為黏性系數。

(5)

(6)

1.2 幾何模型和初邊值條件

1.3 物性參數

本實驗中影響計算結果的物性參數包括線膨脹系數α、Lamé常數λ和μ,以及黏性系數η。其中,線膨脹系數α在整個模型中取為常數10-5K-1[14],而彈性參數和黏性系數則隨半徑變化。彈性參數可以通過下式與月球的密度和波速結構直接關聯：

(7)

其中：Vp和Vs分別為縱波波速和橫波波速。目前,研究者在阿波羅計劃提供的月震數據基礎上使用不同方法計算得到了多種波速結構[19-21],這些波速結構在細節上雖然不盡相同,但反映出的月震波波速隨半徑變化的趨勢是相似的。我們選用Garcia等[20]提出的密度和波速結構(圖1(b)),計算得到的彈性參數隨半徑的變化曲線如圖1(c)所示。

圖1 數值模型的溫度和物性參數分布Fig.1 Distribution of temperature and physical properties of the numerical model

相對于彈性參數,人們對于月球內部介質黏性系數的認知更加模糊。在前人的月球動力學演化模型中,普遍假設月幔的物質組分與地幔相似,滿足Arrhenius定律[14,16-18]

(8)

其中：R為理想氣體常數,Eeff為有效活化能,Tref為參考溫度,η0為Tref下的參考黏度。式(8)對一般情形下的Arrhenius定律做了簡化,省略了環境應力和巖石粒徑對黏度的影響。Zhang等[16]指出,若在式(8)中取Eeff=100～400 kJ/mol,則計算得到的有效黏度與使用標準Arrhenius定律和實驗室條件下測得的橄欖石活化能(335～670 kJ/mol)[22]計算得到的有效黏度相當。對于參考黏度η0,實驗和觀測給出的約束較少,數值實驗中通常取η0=1019～5×1021Pa·s[16-18]。本實驗中,對Eeff=100,200,300 kJ/mol和η0=1019,1020,1021Pa·s的情況分別進行計算,以觀察黏性參數對熱應力的影響。不同流變參數計算得到的黏性系數隨半徑r變化的曲線如圖1(c)所示。

由于月球內部溫度的跨度很大,按式(8)計算會得到不符合實際的黏度值,因此在實際計算中通常將黏度限制在一定范圍內。使用公式

(9)

計算實際黏度,其中黏度上限ηmax在模型中代表月球淺部的黏性系數,直接影響月殼的應力松弛時間。在多數以月幔對流為主的動力學演化模型中,出于數值穩定的需要,通常將ηmax限定在1025Pa·s左右。在Kamata等[11]和Qin等[12]的黏彈性月球模型中,ηmax分別取為1028和1030Pa·s,從而保證月球淺部表現出以彈性為主的力學性質。本實驗中,對ηmax=1026,1028,1030Pa·s的情況分別進行計算,以觀察月球淺部黏性系數對應力松弛的影響。

為方便敘述,下文中使用ExxxMxxCxx的方式描述模型中各個參數的取值,其中E指代有效活化能Eeff,M指代月幔參考黏度η0,C指代黏性系數最大值ηmax(即巖石圈黏度)。例如,模型E100M21C30中各個參數取值為Eeff=100 kJ/mol,η0=1021Pa·s,ηmax=1030Pa·s。同時,省略熱應力的下標thermo,若無特殊說明,下文中的σ均指熱應力。

2 計算結果

2.1 月球巖石圈黏度對熱應力松弛的影響

圖2展示了模型E200M20C26、E200M20C28和E200M20C30的計算結果。由圖2(a)可知,熱應力隨半徑的變化大致分成3段：400 km1 700 km時,σrr和σφφ的絕對值均逐漸減小,在月表處σrr變為0。

圖2 模型E200M20C26、E200M20C28和E200M20C30的計算結果Fig.2 Results of models E200M20C26, E200M20C28, and E200M20C30

圖2的計算結果可以通過簡化的一維模型來解釋：對于一維Maxwell體,當溫度隨時間線性變化時,熱應力一面積累、一面松弛,按下式的形式變化

(10)

其中：t0=η/E為松弛時間,σmax為最終達到平衡時的應力值。式(10)的圖像為：起始時刻應力近似呈線性增長,然后漸漸變緩,最后趨于一個穩定值,即溫度勻速變化時,應力的積累和松弛達到動態平衡。溫度線性變化的速率越大、松弛時間越長,則σmax的絕對值越大,反之σmax的絕對值越小。在我們的數值模型中,月球淺部的黏性系數很高(η≈ηmax),松弛時間很長。當ηmax≥1028Pa·s時熱應力在500 Ma內幾乎一直處在近似線性增加的階段,表現接近彈性體;而月幔深部黏性系數低,松弛時間短,很快達到熱應力積累和松弛接近動態平衡的階段。

(11)

圖2(b)表明,模型E200M20C26的淺部熱應力在500 Ma時基本達到加載和松弛平衡的狀態,而模型E200M20C28和E200M20C30的熱應力仍在勻速加載。實際上,按月幔的剪切模量為70 GPa估算,模型E200M20C26、E200M20C28和E200M20C30的最大松弛時間分別為5 Ma、0.5 Ga和50 Ga,因此后兩者在模型時間尺度內的松弛量可忽略或近于可忽略。Kamata等[11]指出,ηmax≥1028Pa·s的巖石圈對應于干燥斜長石構成的月殼和干燥橄欖石構成的上月幔,而月球的自由空氣重力異常數據支持月球巖石圈具有高黏度。因此,在后面的實驗中統一取ηmax=1030Pa·s。需要注意的是,我們在模型中假定月殼和巖石圈月幔具有相同的黏度,而實際上二者的黏度應當不同。但由于模型淺部的熱應力狀態完全由彈性主導,因此黏度數值1～2個數量級的差別并不影響計算結果。

2.2 活化能和參考黏度對熱應力分布的影響

有效活化能Eeff對熱應力分布的影響如圖3(a)所示。隨著Eeff增大,月幔黏度隨半徑增大的速度加快,導致熱應力持續加載的區域更厚,模型E100M20C30、E200M20C30和E300M20C30中熱應力加載和松弛達到平衡狀態的深度分別為300、450和550 km。同時,更厚的彈性層使內部加載和松弛平衡區域所承受的“靜水壓”更大,模型E100M20C30和E300M20C30的下月幔熱應力相差近2倍。月幔參考黏度η0對熱應力分布的影響如圖3(b)所示。由于ηmax遠大于η0,彈性層厚度幾乎不受η0的影響,因此η0對熱應力的分布影響較小。

圖3 不同模型計算得到的黏性系數、徑向正應力、切向正應力和差應力隨半徑的變化Fig.3 Distribution of viscosity, radial normal stress, tangential normal stress, and differential stress in different models

最后,給出4種端元模型(E100M19C30、E300M21C30、E100M21C30和E300M19C30)的計算結果(圖3(c))。在合理參數范圍內(100 kJ/mol≤Eeff≤300 kJ/mol,1019Pa·s≤η0≤1021Pa·s),熱應力加載和松弛達到平衡狀態的最小和最大深度分別為250和550 km,下月幔所承受的拉張熱應力的最小和最大值分別為10和20 MPa。所有模型均未在下月幔產生超過10 kPa的差應力。

2.3 現今月球內部熱應力的估計

上述計算結果的一個顯著特征是月幔中下部熱應力的加載和松弛相平衡,接近于“靜水壓”狀態,熱應力的大小主要受淺部彈性層的徑向熱應力大小控制,而不受中下部月幔的局部擾動(例如月幔對流)的影響。因此,可以將數值模擬的起始時間前推。

假設月幔發生過重力翻轉事件,并將其視為對熱應力的一次重置。文獻[14,16-18]的計算結果表明月幔翻轉事件的持續時間為0.5～1.5 Ga左右,因此可以假定3 Ga之后不再發生波及整個月幔的對流事件,亦即3 Ga為熱應力積累的起點。使用端元模型E100M19C30和E300M21C30模擬月球3 Ga至今的熱應力演化,計算結果如圖4所示。

圖4 模型E100M19C30和E300M21C30自3 Ga至今的熱應力演化計算結果Fig.4 Thermal stress evolution of models E100M19C30 and E300M21C30 from 3 Ga to present day

由圖4(a)可見,3 Ga時月球內部的溫度分布和溫度變化速率與現今差異明顯。受溫度的影響,3 Ga時月幔的黏性系數遠低于現今(圖4(b)),因此在熱應力演化的早期階段中上月幔的黏性松弛時間很短,導致熱應力加載和松弛達到平衡的深度比500 Ma至今的計算結果更淺(見圖4(c)和4(d),模型E100M19C30和E300M21C30計算得到的熱應力加載和松弛達到平衡的深度分別為200和300 km)。根據一維參量化模型的計算結果,月球巖石圈在3 Ga至今的時間范圍內最薄時約300 km,因此月幔對流不會影響到淺部彈性層的熱應力積累。在無應力釋放的前提下,月表和月殼底部的切向擠壓熱應力分別積累達到150和460 MPa,月幔中下部的拉張熱應力則達到60～90 MPa。值得注意的是,在上述計算中假設月球的彈性參數、有效活化能和參考黏度的分布在3 Ga的時間尺度上保持不變,而在實際的演化過程中月殼和月幔介質的力學性質可能發生更復雜的變化。因此,我們的計算結果可以看作對現今月球內部熱應力分布的估計,進一步的研究依賴于更精細的數值模型和觀測數據。

3 熱應力與月震的聯系

3.1 熱應力與淺源月震的聯系

如引言中所述,月球淺部的切向擠壓熱應力是引發淺源月震的可能因素之一。假設月殼介質遵循Coulomb-Mohr破裂準則,抗剪強度σs可表達為

σs=C+σlithosinθ=C+ρghsinθ,

(12)

其中：C為內聚力,θ為內摩擦角,h表示深度。應力Mohr圓與抗剪強度曲線相切時,有

σdiff/2=ρghsinθ+Ccosθ,

(13)

其中σdiff/2=|σrr-σφφ|/2為應力Mohr圓半徑。為估計熱應力可觸發淺源地震的最大深度,近似取g=1.63 m·s-2,ρ=3 300 kg·m-3,C=1 MPa,θ=30°,此時σdiff/2≥ρghsinθ+Ccosθ的范圍為0～80 km (見圖5黑色虛線與縱軸的交點),即月球淺部0～80 km深度范圍內均有可能發生熱應力導致的壓性剪切破裂。目前對淺源月震的確切深度存在爭議,通過不同波速模型得到的淺源月震最大深度在168～220 km不等[23]。我們的估算結果給出的壓性剪切破裂深度區間大致為淺源地震的震源深度區間的1/3。值得注意的是,內摩擦角θ的取值對壓性剪切破裂的深度區間影響很大,±5°的差別會引起最大破裂深度?10 km左右的變化。

圖5 月球淺部Mohr圓半徑隨深度的變化Fig.5 Variation of the radius of Mohr’s circle against depth

3.2 熱應力與深源月震的聯系

相對于淺源月震,深源月震的發震原因更加難以解釋。盡管深源月震與潮汐作用的關聯得到普遍承認,但潮汐作用引起的差應力非常小(根據張貝等[9]的計算,潮汐作用引起的最大剪應力不超過1 MPa),遠遠達不到3～4 GPa圍壓下巖石的破裂極限。由于波速結構顯示月幔底部可能存在熔融或部分熔融層[24],Frohlich和Nakamura[25]認為處于熔融或部分熔融狀態的月幔物質在潮汐應力的反復作用下發生疲勞,從而可以在較低的差應力下發生破裂。張貝等[9]提出一種不同的思路,認為月幔深部可能存在孔隙結構(孔隙流體即為熔體),巖石骨架承受的有效應力因孔隙壓的存在而減小,從而大大降低了破裂強度。同時,月幔深部的孔隙流體可能在靜巖壓力作用下向上遷移,而在上月幔會受到切向擠壓熱應力的阻礙,從而形成類似于巖漿囊的構造。這一推論與深源月震大多叢集于一些特定區域形成“震源窩”的特征[26]相吻合。

我們的計算結果為深部月?？紫督Y構的發育提供了一種解釋,如圖6所示,圖中plitho表示靜巖壓力,pthermo=-σrr=-σφφ<0表示月球冷卻在月幔深部產生的拉張“靜水壓”,ppore表示孔隙壓,σdiff表示潮汐作用引起的差應力。靜巖壓力的傳遞依靠礦物顆粒的相互擠壓,在微觀尺度上,不同的礦物顆粒由于組分、晶體結構、晶軸方向等等的不同導致力學強度存在差異,因而靜巖壓力的分布可能并不均勻。以圖6(a)所示的情況為例,3顆強度較高的晶粒形成一個結實的“骨架”承受了大部分的靜巖壓力,而位于其中的強度較低的晶粒所承受的靜巖壓力較小。與靜巖壓力不同,熱應力不受礦物顆粒排布的影響,且多數礦物的熱膨脹系數大小相近,因此可以認為熱應力的分布比較均勻。由于深部月幔的拉張熱應力可積累至數十～100 MPa,因此某些局部熱應力的大小可能超過靜巖壓力,從而形成張裂隙。隨著孔隙壓力的增大,巖石的應力Mohr圓向左移(圖6(b)),因此潮汐作用形成的剪應力(～0.5 MPa)也有可能引起月幔深部的剪切破裂,從而觸發深源月震。值得注意的是,靜巖壓力、熱應力和孔隙壓力共同作用下的有效應力可能為拉張狀態,因

圖6 熱應力觸發深源月震的物理機制示意圖Fig.6 A cartoon illustrating the physical mechanism of deep moonquakes triggered by thermal stress

此深源月震的震源機制可能并非走滑,而是張性破裂,這與深源月震的S波和P波振幅分析結果[27]相符。這種解釋隱含的推測是深源月震應力降應該比較小,與實際觀測到的應力降(約為0.05 MPa[8,23])也是吻合的。

4 結論

本研究使用三維黏彈性有限元模型計算月球勻速冷卻的熱應力積累過程。實驗結果表明月球淺部彈性層中的熱應力以切向擠壓為主,而月幔深部的熱應力因黏性松弛而呈現接近“靜水壓”狀態,徑向和切向正應力的大小完全受淺部徑向熱應力的控制,因此可以不受局部月?；顒拥挠绊懚鴱脑虑蜓莼脑缙陔A段積累至今。在高黏度巖石圈假設下,現今月殼底部的切向擠壓應力可積累達到400 MPa,相應的月幔深部拉張應力可達60～120 MPa。我們的計算結果為深源月震震源機制的孔隙流體假設提供了支持,而究竟潮汐應力、熱應力和孔隙流體的共同作用能否使得深部月幔介質達到破裂條件,仍有待于進一步的研究。