崔曉娜 陳林

上地幔部分熔融形成的低粘度、低密度的基性巖漿在地殼中聚集和侵位的過程通常會(huì)伴隨巖漿房、巖墻或巖席的形成(Blacketal., 2021; Caricchietal., 2021)。巖漿的侵位和噴發(fā)往往會(huì)導(dǎo)致地震、海嘯等自然災(zāi)害的發(fā)生(Carvajaletal., 2022; Journeauetal., 2022),也會(huì)形成銅礦、金礦等金屬礦床(Parketal., 2021; 楊林等, 2023)。因此,理解幔源巖漿從深部上升至地表的動(dòng)力學(xué)過程具有重要的科學(xué)意義。高分辨率地球物理成像結(jié)果揭示,我國東北長白山火山下方的地殼內(nèi)部存在多層級(jí)的巖漿聚集區(qū)(Zhuetal., 2019; Yangetal., 2021)。巖石地球化學(xué)數(shù)據(jù)表明,地表出露的玄武巖呈現(xiàn)不同期次的巖漿補(bǔ)給(郭文峰等, 2014; Liuetal., 2021)、巖漿房多次聚集的特征(Leeetal., 2021; Yietal., 2021),且在地表噴發(fā)的火山數(shù)量遠(yuǎn)小于地下賦存巖漿的體積(Crisp, 1984)。由此可見,幔源巖漿在地殼中的侵位具有多期次巖漿供給和分層聚集的特點(diǎn)。然而,幔源巖漿侵位過程中產(chǎn)生的分層聚集樣式受何種因素控制仍不清楚。

針對(duì)幔源巖漿的上涌過程,前人已開展了大量數(shù)值和物理模擬研究。早期數(shù)值模擬研究主要關(guān)注巖漿侵位對(duì)圍巖帶來的熱影響(Annenetal., 2006; Petford and Gallagher, 2001),或是探討單期幔源巖漿在均質(zhì)地殼的聚集過程(Gerya and Burg, 2007; Schubertetal., 2013)。近期的數(shù)值模擬研究主要關(guān)注巖漿的傳輸通道(Gudmundsson, 2006; Rummeletal., 2020)。圍巖的流變強(qiáng)度及其分層性決定了巖漿的傳輸方式,包括熔融底辟、減壓通道或張性裂隙(Kelleretal., 2013),而巖墻和底辟可以共同促進(jìn)巖漿向上傳播(Caoetal., 2016)。伴隨巖漿侵位,侵入體周圍發(fā)育的裂隙由初期呈徑向展布的剪裂隙轉(zhuǎn)換為向上擴(kuò)展的主干張性斷裂(盧靖雯等, 2022),同時(shí)巖漿侵位會(huì)對(duì)圍巖施加附加變形,導(dǎo)致近場與遠(yuǎn)場的小尺度變形(常成和羅綱, 2022)。Gorczyk and Vogt (2018)利用三維熱-力學(xué)模擬方法探究地殼流變和溫度結(jié)構(gòu)對(duì)基性巖侵位形態(tài)的影響,最近Rummeletal.(2020)考慮了巖漿熔融抽取過程中的巖漿成分演化,模擬殼內(nèi)多期次玄武巖席聚集形成巖漿房的機(jī)制。物理模擬研究早期主要探究巖墻向巖席的轉(zhuǎn)換機(jī)制,通過將流體注入固體明膠的模擬實(shí)驗(yàn),驗(yàn)證了當(dāng)固體明膠存在剛度對(duì)比層、層間弱接觸面或外部應(yīng)力時(shí),可以導(dǎo)致垂向傳播的裂隙運(yùn)動(dòng)軌跡轉(zhuǎn)為水平向(Kavanaghetal., 2006,2015; Menandetal., 2010)。最近,Williamsetal.(2022)利用物理模擬方法探究了巖漿沿巖席長距離水平傳播的形成機(jī)制,發(fā)現(xiàn)巖漿的流變學(xué)性質(zhì)具有重要影響。然而,物理模擬受相似比的限制,難以直接運(yùn)用于理解實(shí)際的地質(zhì)過程(Burovetal., 2014)。綜上,前人的相關(guān)研究要么沒有考慮巖漿的多期補(bǔ)給,要么沒有考慮巖墻作為不同層級(jí)巖漿房之間的連接通道,尤其是針對(duì)多期幔源巖漿在地殼中分層侵位過程及其控制因素的研究仍然很薄弱。

在本研究中,我們聚焦幔源巖漿從深部到地殼的聚集運(yùn)移過程,通過在二維熱-力學(xué)數(shù)值模擬程序中引入多期巖漿脈沖和巖墻生成算法,模擬幔源巖漿從深部地幔上升到地殼的動(dòng)力學(xué)過程,系統(tǒng)測試地殼地溫梯度和流變強(qiáng)度的影響,探究幔源巖漿在地殼內(nèi)部分層侵位的控制因素。我們將首先介紹數(shù)值模擬方法與模型設(shè)計(jì);然后描述地殼地溫梯度和地殼強(qiáng)度對(duì)殼內(nèi)巖漿賦存狀態(tài)的影響;最后,討論幔源巖漿在地殼內(nèi)聚集樣式的控制因素。

1 數(shù)值模擬方法

我們使用開源的二維熱-力學(xué)MATLAB程序(Gerya, 2019)模擬幔源巖漿多期次侵位的動(dòng)力學(xué)過程。本程序在求解動(dòng)量守恒方程、質(zhì)量守恒方程與能量守恒方程時(shí),聯(lián)合交錯(cuò)網(wǎng)格有限差分方法和粒子標(biāo)記(marker-in-cell, MIC)技術(shù)。在節(jié)點(diǎn)上解方程,在marker上傳遞物性參數(shù),可以有效避免利用有限差分方法求解連續(xù)性方程時(shí)可能引起的數(shù)值發(fā)散問題。同時(shí),該方法考慮了巖石的粘-彈-塑性流變本構(gòu)關(guān)系和巖石的部分熔融,可以在時(shí)間與空間尺度上精細(xì)模擬巖漿侵位的動(dòng)力學(xué)過程。

1.1 控制方程

(i)質(zhì)量守恒-連續(xù)性方程

(1)

式中,Vx和Vy分別代表速度矢量的水平和垂直分量。

(ii)動(dòng)量守恒-斯托克斯方程

(2)

(3)

(4)

(5)

(iii)能量守恒-熱傳導(dǎo)方程

(6)

1.2 巖石流變

(7)

(8)

(9)

其中,式中,R是氣體常數(shù);AD、Ea、Va和n是巖石物理實(shí)驗(yàn)約束的流變參數(shù),分別代表物質(zhì)常數(shù)、活化能、活化體積和應(yīng)力指數(shù);f為粘度因子,用于線性改變物質(zhì)的有效粘度(Beaumontetal., 2006)。由于實(shí)驗(yàn)室條件下物質(zhì)成分的變化,所測得的流變參數(shù)也具有很大的不確定性。為了將這些結(jié)果應(yīng)用到實(shí)際地質(zhì)過程,我們使用粘度因子對(duì)物質(zhì)在實(shí)驗(yàn)室條件下測得的有效粘度進(jìn)行線性調(diào)節(jié)(劉丹紅和陳林,2020)。在這里,我們通過改變粘度因子調(diào)控地殼的強(qiáng)度。

熔融巖石(5%

表1 數(shù)值模型中用到的主要物理參數(shù)

物質(zhì)的塑性變形遵從Drucker-Prager屈服準(zhǔn)則(Ranalli,1995):

σyield=C0+Psin(φ)

(10)

其中C0是巖石的內(nèi)聚力;φ是有效內(nèi)摩擦角。

1.3 部分熔融

根據(jù)巖石物理實(shí)驗(yàn)數(shù)據(jù)的約束,程序考慮了巖石的部分熔融行為(Gerya and Yuen, 2003a; Burg and Gerya, 2005)。假設(shè)部分熔融百分?jǐn)?shù)M與溫度為簡單的線性關(guān)系:

(11)

式中,Tsolidus為巖石的固相線,Tliquidus為巖石的液相線,各種巖石的固相線和液相線見表1。

部分熔融巖石的有效密度ρeff取決于部分熔融百分?jǐn)?shù)M:

ρeff=ρsolid-M(ρsolid-ρmolten)

(12)

式中,ρsolid和ρmolten分別代表巖石的固相密度和熔融密度,它們都是關(guān)于溫度和壓力的函數(shù),計(jì)算表達(dá)式見表1。其中,我們對(duì)幔源巖漿密度的選取參照了Sanloup (2016)的結(jié)果,并且與該作者匯編的結(jié)果進(jìn)行了對(duì)比,在0～5GPa的壓力范圍內(nèi),我們的結(jié)果中巖漿的密度與前人的實(shí)驗(yàn)結(jié)果和理論計(jì)算結(jié)果均一致(電子版附圖1),說明巖漿密度的選取具有合理性。

圖1 巖墻形成的示意圖

部分熔融對(duì)有效熱容量(Cpe)和熱擴(kuò)散系數(shù)(αe)的影響由下式給出(Burg and Gerya, 2005):

(13)

1.4 多期巖漿脈沖的實(shí)現(xiàn)

火山噴發(fā)需要體積足夠大、低粘度、低密度、部分熔融且持續(xù)供應(yīng)的巖漿熱源,以提供浮力、維持高溫(Caricchietal., 2021)。我們通過預(yù)設(shè)巖石圈之下的巖漿房來模擬幔源巖漿上升的動(dòng)力學(xué)過程,在已有的二維粘-彈-塑性數(shù)值模擬程序(Gerya and Burg, 2007)的基礎(chǔ)上加入多期巖漿脈沖算法,在巖石圈以下多期次、定量的穩(wěn)定提供巖漿熱源,盡量保證巖漿通道不變形。算法實(shí)現(xiàn)細(xì)節(jié)為:(1)我們添加新一期巖漿脈沖的判斷條件為:深度120km以下,部分巖漿通道和巖漿源區(qū)中熔體的總含量小于5%;(2)當(dāng)熔體含量滿足上述條件時(shí),我們用性質(zhì)相同但標(biāo)號(hào)不同的marker補(bǔ)充新一期次的幔源巖漿;(3)在模擬過程中,我們總共提供六期巖漿脈沖,為巖漿的持續(xù)上涌提供動(dòng)力。

1.5 巖墻的生成算法

巖墻作為拉伸作用形成的薄弱帶,它的形成與圍巖的應(yīng)力場有關(guān)(Gudmundsson, 2006; Bertelsenetal., 2018)。在本文中,拉伸應(yīng)力標(biāo)記為正、擠壓應(yīng)力標(biāo)記為負(fù)。巖墻作為巖漿房之間的連接通道(Cox, 1980; Reuberetal., 2018),巖漿沿巖墻傳播(Becerriletal., 2013; Scuderoetal., 2019),其寬度隨深度變化(Geshietal., 2010, 2020)。目前對(duì)于巖墻的熱-力學(xué)模擬主要探究巖墻對(duì)巖漿上升過程的影響,例如,Caoetal.(2016)認(rèn)為花崗巖巖漿房會(huì)形成輻射狀巖墻,為了探究巖漿上升過程中巖墻和底辟的相互作用,因此設(shè)定臨界的應(yīng)力和第二應(yīng)力不變量作為垂直狀巖墻形成的條件。Rummeletal.(2020)通過巖石熱-力學(xué)模擬方法,考慮巖墻作為玄武巖熔融提取的通道。我們通過探究幔源巖漿侵位對(duì)地殼應(yīng)力狀態(tài)帶來的影響,認(rèn)為垂直狀巖墻是巖漿上升的通道,輻射狀巖墻是巖漿房體積擴(kuò)張時(shí)圍巖的同期破裂,表征巖漿房的發(fā)育。因此,本文中考慮巖墻是在巖漿侵位對(duì)上覆圍巖造成拉伸環(huán)境中形成的(Gudmundsson, 1988),作為巖漿向淺部運(yùn)移的通道。

我們?cè)诔绦蛑屑尤霂r墻生成算法,巖墻的形成和擴(kuò)展由巖漿房超壓值和地殼應(yīng)力共同控制。巖漿房超壓值為巖漿房內(nèi)部的總壓強(qiáng)與靜巖壓力的差值(巖漿房內(nèi)部的超壓值相差無幾),巖漿房內(nèi)部的超壓值會(huì)隨著每一期次的巖漿補(bǔ)給而增加、隨著巖漿房的體積增大而減小,超壓值演化表征了巖漿房的發(fā)育過程;當(dāng)巖漿底侵時(shí),巖漿聚集使得巖漿房之上產(chǎn)生應(yīng)力集中的區(qū)域,即巖漿房的發(fā)育對(duì)上覆的上地殼與下地殼形成了水平的拉伸應(yīng)力,這種拉伸應(yīng)力即為地殼應(yīng)力,地殼的應(yīng)力場決定了巖墻垂直向發(fā)育。以下地殼巖墻的形成過程為例,具體算法如下:(1)首先,由于多期次的巖漿脈沖供應(yīng),殼內(nèi)巖漿房發(fā)育至一定體積,將巖漿房內(nèi)部超壓值小于120MPa作為巖漿房發(fā)育結(jié)束的時(shí)刻與判斷是否形成巖墻的初始時(shí)刻。(2)隨后,當(dāng)巖漿房之上的地殼應(yīng)力大于30MPa時(shí),發(fā)育長度為巖漿房頂部至上、下地殼界面、寬度為2km的巖墻(見圖1)。經(jīng)過多次測試,2km寬度的巖墻是保證巖漿上升的最窄寬度,過窄的巖墻會(huì)導(dǎo)致通道快速冷卻、無法輸運(yùn)巖漿,因此我們?cè)谀M中設(shè)定巖墻的最小寬度為2km。

2 模型設(shè)計(jì)

為了模擬幔源巖漿上升的動(dòng)力學(xué)過程,我們?cè)O(shè)計(jì)了如圖2所示的初始模型。模型在水平方向上的寬度為160km,垂向深度為180km,包含20km上地殼、20km下地殼和90km巖石圈地幔,代表了一個(gè)典型的大陸巖石圈,具體物性參數(shù)見表1。我們初始在巖石圈底部預(yù)設(shè)一個(gè)半徑為10km、溫度為1700K的圓形幔源巖漿區(qū)域,在頂部設(shè)置了一個(gè)寬度為3km的巖漿通道,連接到地殼底部。整個(gè)模型空間被321×361個(gè)節(jié)點(diǎn)組成的規(guī)則網(wǎng)格離散化,空間分辨率為500m,總共120萬個(gè)粒子(markers)隨機(jī)分布在網(wǎng)格空間中,粒子分辨率高達(dá)150m。

圖2 模型設(shè)計(jì)

模型的所有力學(xué)邊界均為自由滑邊界。在模型的上部,我們?cè)O(shè)計(jì)一個(gè)厚20km的偽空氣層,其密度和粘度分別為1kg/m3、1013Pa·s,用于模擬巖漿上升過程中產(chǎn)生的地表起伏(Gerya and Burg, 2007; Schmelingetal., 2008)。模型的熱邊界條件:頂邊界固定為273K,側(cè)邊界為絕熱邊界,底邊界為固定為1600K (Burg and Gerya, 2005; Huetal., 2018),整個(gè)空氣層溫度固定為273K,巖石圈底邊界以下為絕熱地溫梯度0.5K/km。為了區(qū)分地殼的不同熱狀態(tài),我們采用分段函數(shù)賦予了地殼兩種初始溫度結(jié)構(gòu):第一種溫度結(jié)構(gòu)初始固定莫霍面的溫度,設(shè)定地表和莫霍底界面的初始溫度分別為273K和TMoho=900K,通過改變上、下地殼界面的初始溫度(TUC=586K、686K、786K)來調(diào)節(jié)上、下地殼的地溫梯度(GUC、GLC),讓溫度隨深度線性增加;第二種溫度結(jié)構(gòu)初始固定巖石圈底部的溫度,設(shè)定地表和巖石圈底界面的初始溫度分別為273K和1573K,通過改變莫霍面的初始溫度(TMoho=673K、773K、873K、973K、1073K)來調(diào)節(jié)地殼和巖石圈的地溫梯度,同樣讓溫度隨深度線性增加。

3 模擬結(jié)果

為了探究來自地幔深部的多期巖漿在地殼中分層聚集的動(dòng)力學(xué)過程,我們基于多期巖漿脈沖和巖墻算法的MATLAB程序,開展了27組數(shù)值模擬實(shí)驗(yàn)(表2)。在模擬中,系統(tǒng)測試了地殼地溫梯度和地殼強(qiáng)度對(duì)巖漿上涌過程的影響。在這里,我們選擇上地殼的粘度因子fUC=0.05、下地殼的粘度因子fLC=0.2、TUC=686K、TMoho=900K的模型作為參考模型(即表2中的model15)。

表2 數(shù)值模擬實(shí)驗(yàn)

3.1 參考模型

圖3為參考模型(model15,fUC=0.05,fLC=0.2,GUC=20K/km,GLC=10K/km)的演化結(jié)果。來自巖石圈底部的幔源巖漿穿過巖石圈到達(dá)莫霍面,初始階段(4.8kyr;圖3ai),巖漿在莫霍面之上聚集形成橢圓狀巖漿房,伴隨同時(shí)的垂向增厚和橫向擴(kuò)展,粘度結(jié)果顯示巖漿房周圍出現(xiàn)局部的低粘度區(qū)域(圖3aii),表征巖漿房的擴(kuò)張趨勢。隨后,由于巖漿侵位對(duì)上覆地殼造成的拉伸力,下地殼形成巖漿通道使得巖漿由此上升(圖3bi),水平拉伸應(yīng)力瞬時(shí)消散(圖3biii)。8.6kyr時(shí),巖漿上涌至上、下地殼界面,多期次巖漿脈沖帶來的高超壓導(dǎo)致巖漿通道不斷變寬,在巖漿房內(nèi)部,二期巖漿及多期巖漿以韌性底辟形式對(duì)流, 促進(jìn)多期巖漿混合及巖漿房上升(圖3ci),粘度結(jié)果顯示巖漿房之上局部的低粘度區(qū)域(圖3cii),表征巖漿的上升趨勢。最終13.1kyr時(shí),上、下地殼界面之上巖漿聚集形成扁平狀巖漿房,但巖漿房頂部的應(yīng)力累積不足(圖3ciii),達(dá)到穩(wěn)態(tài)至形成雙層聚集的巖漿房,此時(shí)下地殼巖漿房體積擴(kuò)張,并且內(nèi)部卷入極少量的地殼和巖墻等物質(zhì)(圖3di),上地殼與下地殼減薄程度不同的原因在于兩者對(duì)于巖漿侵位帶來的空間問題響應(yīng)不同。在13.1kyr的時(shí)間尺度內(nèi),巖漿房內(nèi)部保持活躍,進(jìn)行持續(xù)的巖漿對(duì)流。

圖3 參考模型

圖3中aiii～diii為參考模型的水平偏應(yīng)力隨時(shí)間的變化,下地殼與上地殼巖漿侵位均對(duì)巖漿房之上的局部區(qū)域形成水平拉伸應(yīng)力。上地殼累積的拉伸應(yīng)力在脆韌性轉(zhuǎn)換深度處最大,向地表逐漸降低,向深部為韌性區(qū)應(yīng)力累積很弱;下地殼累積的拉伸應(yīng)力在上、下地殼界面處最大,向深部減小,同樣韌性區(qū)應(yīng)力累積很弱,并且下地殼巖漿侵位導(dǎo)致上、下地殼界面隆起,可以被上地殼韌性區(qū)容納(圖3bii)。同時(shí),我們選取參考模型(model15)中三個(gè)固定點(diǎn)的水平偏應(yīng)力結(jié)果來描述巖漿上涌對(duì)下地殼、上地殼、地表帶來的影響,在圖3aiii中用紅色的圓標(biāo)記。下地殼巖漿房之上紅色固定點(diǎn)拉伸力的初始波峰表示下地殼巖漿聚集的初始階段,第二個(gè)40MPa小波峰表示4.5kyr添加的第二期巖漿脈沖,4.8kyr時(shí)拉伸應(yīng)力由38MPa驟降表征巖漿房破裂,0MPa表征該深度已經(jīng)轉(zhuǎn)變?yōu)閹r墻或巖漿成分(見圖1)。與之相似,上地殼脆韌性轉(zhuǎn)換深度處,應(yīng)力的第一個(gè)波峰表示下地殼巖漿聚集過程中上地殼的同期應(yīng)力累積,第二個(gè)34MPa小波峰表示添加的第二期巖漿脈沖,第三個(gè)41MPa的波峰表示上地殼巖漿的聚集過程,隨后應(yīng)力逐漸消散。巖漿初始聚集時(shí)地表應(yīng)力的累積滯后,隨后在巖漿開始沿巖墻向上傳播至源區(qū)六期巖漿脈沖結(jié)束的過程中,淺表一直處于拉伸環(huán)境,當(dāng)巖漿脈沖停止后,應(yīng)力出現(xiàn)負(fù)值表示地殼淺部轉(zhuǎn)為擠壓狀態(tài)。

我們選擇下地殼巖漿房底部與上地殼巖漿房底部的兩個(gè)固定點(diǎn)從超壓演化分析巖漿房的發(fā)育過程,在圖3di中用白色的圓標(biāo)記。下地殼巖漿房內(nèi)的超壓因一期巖漿初始聚集而驟增(圖4a),后因巖漿房的不斷發(fā)育,超壓值隨著巖漿房體積增大而持續(xù)減小,后五期的巖漿脈沖均可引起超壓值驟增,使巖漿房維持高超壓狀態(tài)而體積不斷擴(kuò)張,4.8kyr時(shí)超壓值驟降表征巖漿房的破裂(見圖1),此時(shí)為下地殼巖漿房發(fā)育結(jié)束的時(shí)刻與形成下地殼巖墻的時(shí)刻,最終下地殼巖漿房內(nèi)部的超壓值趨于穩(wěn)定,當(dāng)巖漿脈沖停止后超壓負(fù)值表征巖漿上升停滯。上地殼巖漿房內(nèi)部的超壓初始值對(duì)應(yīng)下地殼巖漿房破裂的時(shí)刻(圖4b),超壓值由負(fù)轉(zhuǎn)正表征上、下地殼界面深度的初始巖漿聚集,后四期巖漿脈沖引起的超壓最大值小于30MPa,當(dāng)巖漿脈沖停止后超壓負(fù)值亦表征巖漿上升停滯。

圖4 巖漿房內(nèi)超壓值隨時(shí)間的演化

3.2 地殼地溫梯度的影響

在參考模型的基礎(chǔ)上,我們測試了地殼地溫梯度的影響。如圖5d,設(shè)定莫霍面的初始溫度分別為TMoho=783K、873K、973K (model23、model24、model25;見表2),根據(jù)單一變量測試原則,其余參數(shù)設(shè)置與參考模型相同。首先,我們對(duì)比13.1kyr時(shí)冷地殼、溫地殼與熱地殼背景下巖漿的賦存狀態(tài)。圖5ai中,冷地殼的地溫梯度為GUC=GLC=12.5K/km,低的地溫梯度使得地殼的粘度過大(5aii),即使巖漿高超壓也無法克服地殼的強(qiáng)流變實(shí)現(xiàn)巖漿聚集,此時(shí)巖漿主要在巖石圈底部聚集。圖5bi中,溫地殼的地溫梯度為GUC=GLC=15K/km,在地殼具有較高的地溫梯度時(shí),巖漿在莫霍面之上聚集形成橢圓狀巖漿房,并且?guī)r漿可以上升至上、下地殼界面。圖5ci中,熱地殼的地溫梯度為GUC=GLC=17.5K/km,高的地溫梯度使得地殼的粘度降低,此時(shí)巖漿僅在莫霍面之上聚集形成席狀巖漿房,產(chǎn)生以上現(xiàn)象是基于巖漿的正浮力與相應(yīng)地殼地溫梯度下的地殼強(qiáng)度共同作用的結(jié)果。

圖5 地殼地溫梯度對(duì)巖漿侵位深度的影響

我們進(jìn)一步分析了巖漿侵位對(duì)地形帶來的影響,圖6a-c分別為冷地殼、溫地殼、熱地殼背景下地形隨時(shí)間的變化,溫地殼與熱地殼的地形演化呈現(xiàn)相同的趨勢:侵位均導(dǎo)致地表抬升,初始中心隆起后轉(zhuǎn)變?yōu)橹行陌枷輧蓚?cè)隆起形態(tài)(圖6b, c),而冷地殼時(shí),地表的地形無明顯起伏(圖6a)。速度驟增均對(duì)應(yīng)多期巖漿脈沖的巖漿供應(yīng)時(shí)刻(圖6f)。圖6e為 13.1kyr時(shí)的地形垂直剖面,地表最大隆起值發(fā)生在溫地殼下,地表最高高程約為3km,熱地殼次之,地表最高高程約為1.5km,冷地殼時(shí),地表的相對(duì)高差小于140m。

圖6 地殼地溫梯度對(duì)地形的影響

綜上,地殼地溫梯度對(duì)巖漿的侵位深度具有重大影響,溫地殼有利于巖漿上升到更淺的深度,地殼的地溫梯度過高或過低都不利于巖漿的上升。并且,在巖漿多期次侵位的過程中,地形多期次抬升,殼內(nèi)巖漿房的高度對(duì)地表的隆起高度貢獻(xiàn)更大。

3.3 地殼強(qiáng)度的影響

為了探究不同的地殼強(qiáng)度對(duì)巖漿聚集的影響,我們引入了粘度因子方法調(diào)節(jié)地殼強(qiáng)度,在參考模型的基礎(chǔ)上,系統(tǒng)測試了fUC=1、0.5、0.1、0.05、0.01,fLC=1、0.5、0.2、0.1時(shí)的地殼強(qiáng)度對(duì)巖漿上升過程的影響(model2～21;見表2)。在這里,我們展示model2、model17、model6三個(gè)端元模型的演化結(jié)果,其他模擬結(jié)果見電子版附圖2。圖7ai為fUC=1、fLC=1(model2)的模擬結(jié)果,在地殼強(qiáng)流變的情況下,巖漿主體在巖石圈底部聚集,僅有小體積的巖漿在下地殼底部侵位。當(dāng)固定上地殼強(qiáng)度、僅變?nèi)跸碌貧?qiáng)度時(shí)(fUC=1,fLC=0.1,model17,圖7bi),巖漿在莫霍面之上形成了扁平狀巖漿房,外部為老期次巖漿內(nèi)部為新期次巖漿,上、下地殼界面發(fā)生寬緩的隆升。當(dāng)固定下地殼強(qiáng)度、僅變?nèi)跎系貧?qiáng)度時(shí)(fUC=0.01,fLC=1,model6,圖7ci),巖漿傾向于向上傳播,破裂前的下地殼發(fā)育球狀巖漿房,巖漿由巖墻穿過下地殼后(見圖1),高超壓驅(qū)動(dòng)下的巖漿在上地殼自發(fā)上升至地表,上、下地殼界面處的韌性區(qū)僅有少量巖漿聚集,地表噴發(fā)主要為新期次的巖漿,包裹極少量的下地殼巖墻、上地殼巖石等物質(zhì)。

圖7 地殼強(qiáng)度對(duì)巖漿上升的影響

綜上,我們基于不同的上、下地殼強(qiáng)度背景下模擬巖漿上涌的過程,表明上、下地殼的相對(duì)強(qiáng)度對(duì)巖漿的聚集樣式具有重大影響,地殼強(qiáng)度對(duì)巖漿房的橫向?qū)挾扔绊懜?地殼越弱,巖漿房的寬度越寬。

4 討論

4.1 巖漿聚集樣式的控制因素

根據(jù)模型測試,我們發(fā)現(xiàn)有兩種因素會(huì)影響殼內(nèi)巖漿的聚集樣式。首先是地殼的地溫梯度,地殼越冷,巖漿無法在殼內(nèi)聚集,地殼越熱,巖漿越易在莫霍面之上聚集;其次是上、下地殼的相對(duì)強(qiáng)度,巖漿主體在流變較弱的深度聚集。在模擬過程中,我們?cè)O(shè)定來自地幔深部的巖漿熱源提供六期巖漿脈沖,添加到巖石圈內(nèi)部的巖漿體積是相同的,排除了由于巖漿體積不同帶來的影響,因此巖漿的聚集樣式主要取決于地殼的強(qiáng)度和熱狀態(tài)。

對(duì)于巖漿聚集樣式的控制因素,Gerya and Burg (2007)通過模擬單期幔源巖漿在均質(zhì)下地殼聚集的過程,強(qiáng)調(diào)圍巖內(nèi)聚力和內(nèi)摩擦系數(shù)在控制侵入體的高度和形態(tài)上起著重要作用,韌性的下地殼有利于巖漿的橫向擴(kuò)張,塑性下地殼使得巖漿沿?cái)鄬由嫌恐恋乇?。Schubertetal.(2013)通過模擬鎂鐵質(zhì)巖漿注入預(yù)設(shè)長英質(zhì)巖漿房的過程,發(fā)現(xiàn)地殼的塑性強(qiáng)度決定了巖漿上升的高度,巖漿侵位低強(qiáng)度的地殼僅會(huì)停留在預(yù)設(shè)的下地殼巖漿房,巖漿侵位中等強(qiáng)度的地殼可以沿?cái)鄬由仙辽系貧?巖漿侵位高強(qiáng)度的地殼僅會(huì)上升至上、下地殼界面。Gorczyk and Vogt (2018)利用三維熱-力學(xué)模擬的方法探究地殼溫度和流變對(duì)巖漿侵位形態(tài)的影響,發(fā)現(xiàn)地殼的溫度越高,巖漿更易在莫霍面聚集形成巖漿房,地殼的塑性強(qiáng)度越小,巖漿侵位的深度越淺。但他們的結(jié)果均沒有形成分層聚集的巖漿房。我們的模擬結(jié)果與前人的研究結(jié)果基本一致,但我們通過引入粘度因子來調(diào)節(jié)地殼的粘性流變性質(zhì),發(fā)現(xiàn)地殼粘性流變分層性在控制巖漿多層級(jí)的聚集過程中具有重要影響。

以參考模型的地溫梯度為基礎(chǔ),我們系統(tǒng)地測試了不同上地殼、下地殼強(qiáng)度時(shí)巖漿的上升過程,并將模擬結(jié)果總結(jié)在圖8中。結(jié)果顯示,當(dāng)上地殼與下地殼均具有較強(qiáng)的流變時(shí)(fUC≥1、fLC≥1),幔源巖漿主體在巖石圈底部聚集;當(dāng)下地殼的強(qiáng)度較弱(fUC<1、fLC<1)時(shí),在莫霍界面處可以形成下地殼巖漿房;進(jìn)一步,當(dāng)上地殼的強(qiáng)度較弱時(shí)(fUC≤0.05、fLC≤0.5),在上、下地殼界面和莫霍面可以形成雙層聚集的上地殼巖漿房和下地殼巖漿房。巖漿在巖石圈、莫霍面、上、下地殼界面等不同深度聚集與上、下地殼的強(qiáng)度呈現(xiàn)了一定的趨勢,強(qiáng)度越弱巖漿房越扁平,上地殼巖漿房更難形成,并且上地殼巖漿房的厚度總是小于下地殼巖漿房的厚度。

圖8 上、下地殼強(qiáng)度對(duì)巖漿聚集樣式的影響在這些模式中,上地殼與下地殼的粘度因子是不同的,觀測到巖漿聚集的三種模式:僅在巖石圈聚集、僅在下地殼聚集與雙層聚集

4.2 對(duì)長白山火山活動(dòng)的啟示

長白山火山歷史上曾發(fā)生多次劇烈噴發(fā)(Weietal., 2013),目前該火山仍存在著較高的潛在噴發(fā)危險(xiǎn)(劉嘉麒等, 2015; Zhangetal., 2018)。前人運(yùn)用了不同的地震學(xué)探測方法探究了長白山火山區(qū)下方的地殼結(jié)構(gòu),結(jié)果顯示,該火山下方的地殼內(nèi)存在多層聚集的巖漿房,在下地殼表現(xiàn)為球狀巖漿房、在上地殼表現(xiàn)為扁平狀巖漿房(張先康等, 2002; Songetal., 2007; Choietal., 2013; 陳棋福等,2019; Fanetal., 2022)。大地電磁結(jié)果顯示該地區(qū)火山下的地殼存在顯著的電導(dǎo)率不均一的特征,對(duì)應(yīng)為巖石成分或熔流體含量的不均一(湯吉等, 2001; 仇根根等, 2014; Yangetal., 2021),預(yù)示著殼內(nèi)存在多層級(jí)的流變性質(zhì)分界面。來自地幔深部的低波速、高溫物質(zhì)為長白山地區(qū)的巖漿活動(dòng)持續(xù)提供熱源(Lei and Zhao, 2005; Tangetal., 2014; Lietal., 2020),維持了長白山火山現(xiàn)今較熱的地溫狀態(tài)(平均熱流值大于70mW/m2,Jiangetal., 2019)。我們的模擬結(jié)果表明,當(dāng)?shù)貧鷰r流變分層性越顯著、地溫梯度越高時(shí),更容易形成多層級(jí)聚集的巖漿房,長白山殼內(nèi)結(jié)構(gòu)特征與附圖2biv所示的結(jié)果具有很好的一致性。因此,我們推測長白山地區(qū)殼內(nèi)巖漿的聚集樣式主要受圍巖地溫梯度和地殼內(nèi)強(qiáng)烈的流變分層性的影響。

4.3 幔源巖漿侵位過程中地殼熱狀態(tài)對(duì)巖漿活動(dòng)區(qū)地形演化的影響

在這里,我們通過設(shè)定不同的背景地殼地溫梯度模擬了幔源巖漿的上升過程,進(jìn)一步分析地表的地形演化結(jié)果發(fā)現(xiàn):在不同的背景地殼地溫梯度下,巖漿活動(dòng)區(qū)的地形演化具有明顯的分類特征,冷地殼背景下巖漿底侵對(duì)殼內(nèi)及地表基本沒有影響,而溫地殼和熱地殼背景下,由于殼內(nèi)巖漿多期次聚集形成巖漿房,巖漿房之上的地表地形出現(xiàn)了多期次的抬升。這種地殼熱狀態(tài)與地表高程之間的對(duì)應(yīng)關(guān)系與許多地質(zhì)實(shí)例相吻合。例如,我國大同火山的背景熱流值約為50～60mW/m2(Jiangetal., 2019),屬于中低熱流區(qū),對(duì)應(yīng)為我們的模擬中的冷地殼背景,該地區(qū)的地形相對(duì)高差較小(電子版附圖3a),與我們針對(duì)冷地殼的地形模擬結(jié)果具有高度的一致性。長白山火山區(qū)的背景熱流值約為70～80mW/m2(Lucazeau, 2019;附圖3b),對(duì)應(yīng)為我們的模擬中的溫地殼背景,相對(duì)于周圍的地形,該火山區(qū)地形的相對(duì)高差約2km(附圖2);黃石超級(jí)火山是美國西部的一個(gè)地幔柱成因火山(Huangetal., 2015),該地區(qū)的平均熱流值高于100mW/m2(Lucazeau, 2019),為典型的熱地殼,黃石超級(jí)火山區(qū)作為巖漿底侵區(qū)域(Colónetal., 2018),相對(duì)于斯內(nèi)克河平原,相對(duì)高差約1.5km(Peng and Humphreys, 1998, 附圖3c)。這兩個(gè)火山區(qū)均由于巖漿底侵導(dǎo)致了地表抬升,與我們針對(duì)溫地殼和熱地殼的地形模擬結(jié)果高度吻合。因此,在幔源巖漿侵位的過程中,地殼的熱狀態(tài)與地表的地形演化之間具有極大的相關(guān)關(guān)系,結(jié)合我們對(duì)不同地殼地溫梯度背景下巖漿侵位過程中巖漿活動(dòng)區(qū)地形的演化結(jié)果,可以得出:巖漿活動(dòng)區(qū)的地形演化能夠反映巖漿侵位時(shí)地殼的熱狀態(tài),地殼越冷,地形的相對(duì)變化越小,地殼越熱,殼內(nèi)聚集的巖漿房會(huì)在地表形成相應(yīng)的地形抬升,抬升幅度與巖漿房的高度正相關(guān)。

4.4 模型局限性

我們采用高分辨率的二維熱-力學(xué)模型模擬了幔源巖漿上升的精細(xì)動(dòng)力學(xué)過程。但是,我們的模型也存在一些局限性,首先是巖漿的粘度,由于巖漿熔流體含量的復(fù)雜性,我們考慮部分熔融巖石的流變時(shí)對(duì)粘度做了一定的簡化,這與實(shí)際地質(zhì)過程中熔融巖石的粘度范圍存在差別(Dingwell, 2006),熔融巖石的粘度與溫度和壓力有關(guān),模擬過程中綜合考慮以上因素會(huì)對(duì)模擬結(jié)果帶來影響,但是不會(huì)影響我們對(duì)殼內(nèi)巖漿的聚集樣式的探究。其次是巖墻的寬度,實(shí)際地質(zhì)過程中巖墻的寬度為厘米或者米級(jí)尺度(Geshietal., 2010, 2020),受計(jì)算條件的限制,我們模擬的巖墻寬度超過了實(shí)際的分辨率。雖然我們的模型不能解釋自然界所有的復(fù)雜性,但是我們可以在極短時(shí)間內(nèi)模擬更為接近實(shí)際的過程,并給定了巖漿上升的一般性模型。未來我們將添加兩相流的模擬(Kelleretal., 2013; R?ssetal., 2019),考慮揮發(fā)分出溶帶來的影響(Tian and Buck, 2022)。

5 結(jié)論

我們通過二維熱-力學(xué)方法模擬了幔源巖漿上涌的動(dòng)力學(xué)過程,探究了地殼地溫梯度和地殼強(qiáng)度對(duì)巖漿上升及其聚集過程的影響,主要獲得以下結(jié)論:

(1)地殼地溫梯度對(duì)巖漿的侵位深度有重要影響。巖漿侵位冷地殼,巖漿主體在巖石圈深度聚集;巖漿侵位溫地殼形成下地殼巖漿房,并且?guī)r漿可以侵位到上、下地殼界面;巖漿侵位熱地殼僅形成下地殼席狀巖漿房。同時(shí),伴隨巖漿脈沖的多期次侵位,地表呈現(xiàn)多期次快速抬升的特征。在溫地殼與熱地殼條件下,巖漿房上方的地表地形呈現(xiàn)中心凹陷兩側(cè)隆起的形態(tài);在冷地殼條件下,地表基本無起伏。

(2)地殼流變性質(zhì)的分層性決定了巖漿侵位的樣式。同一地殼地溫梯度下,強(qiáng)的上地殼與強(qiáng)的下地殼(fUC≥1、fLC≥1)使得巖漿主體在巖石圈聚集,無法上升至地殼;當(dāng)上地殼的強(qiáng)度強(qiáng)、下地殼的強(qiáng)度弱時(shí)(fUC≥1、fLC<1),巖漿在下地殼聚集形成巖漿房;當(dāng)上地殼的強(qiáng)度弱、下地殼的強(qiáng)度強(qiáng)時(shí)(fUC≤0.05、fLC≥1),巖漿傾向于在地表噴發(fā);在較弱的上地殼與較弱的下地殼的情況下(fUC≤0.05、fLC≤0.5),巖漿在殼內(nèi)更容易形成雙層聚集的巖漿房。此外,我們還發(fā)現(xiàn)在殼內(nèi)巖漿多層聚集過程中,周期性巖漿脈沖的補(bǔ)給會(huì)導(dǎo)致巖漿房內(nèi)部超壓值呈現(xiàn)瞬時(shí)驟增隨后下降的特征。

(3)根據(jù)我們的分析,長白山火山殼內(nèi)巖漿的分層聚集樣式主要受控于地殼的熱狀態(tài)和流變性質(zhì)分層性;火山區(qū)地形的相對(duì)變化能夠反映巖漿侵位時(shí)地殼的熱狀態(tài),地殼越冷,地形的相對(duì)變化越小,地殼越熱,殼內(nèi)聚集的巖漿房會(huì)在地表形成顯著的地形抬升。我們的模擬結(jié)果對(duì)理解巖漿在殼內(nèi)的聚集樣式和巖漿活動(dòng)區(qū)的地形演化具有啟示意義。

致謝本研究使用的MATLAB源程序來自www.cambridge.org/9781107143142,所有的圖片均使用GMT繪制(https://www.soest.hawaii.edu/gmt/; Wesseletal., 2019)。感謝三位審稿專家和編輯提出的建設(shè)性意見,這些意見極大地幫助我們提高了文章的質(zhì)量。感謝北京大學(xué)的李揚(yáng)教授、中國科學(xué)院地質(zhì)與地球物理研究所的劉博達(dá)副研究員對(duì)本文前期的研究工作提供的建議,感謝課題組顏智勇、唐嘉萱、謝仁先、司潤港、向宵對(duì)本文的幫助。感謝國家超級(jí)計(jì)算天津中心提供“天河一號(hào)”計(jì)算平臺(tái),感謝北京超級(jí)云計(jì)算中心(BSCC)提供的高性能計(jì)算資源。