地殼中巖漿沿已有斷層運移的數(shù)值模擬——在長白山天池火山的應(yīng)用

程旭，陳祖安，白武明

1中國科學院地質(zhì)與地球物理研究所地球深部研究重點實驗室，北京 100029

2中國科學院大學，北京 100049

1 引言

幔源的巖漿在上升侵入地殼之后，會在不同深度形成巖漿囊或巖漿房.這些巖漿囊與下部的深層巖漿房仍有通道連接.當不斷有深部巖漿補給到淺部的巖漿囊中時，巖漿囊中的壓力就會持續(xù)增大，由此可能會導致巖漿運移.巖漿在地殼中的運移分為兩種類型：一種是巖漿壓力驅(qū)動下圍巖不斷破裂，巖脈不斷增長；另一種是巖漿沿著已有斷層或通道運移.

巖漿運移機制及與之相關(guān)的圍巖破裂是一個巖漿從開始緩慢侵入逐漸過渡到急速上涌噴發(fā)的過程.要全面認識這一過程，除了需要先進的儀器來監(jiān)測深部應(yīng)力外，解析分析與數(shù)值模擬流體與固體破裂相耦合的這一動力學過程也是相當重要的.關(guān)于巖漿在地殼中運移過程的力學已經(jīng)取得了許多研究成果.關(guān)于從巖漿囊開始的巖脈傳播的研究大部分集中在線彈性斷裂力學框架上的斷裂擴展過程，這些同時也考慮了巖漿流動的流體動力學過程.Weertman（1971）考慮了一個受到張應(yīng)力和重力誘發(fā)的靜水壓力下的水平彈性板中含流體裂紋的問題，研究結(jié)果表明足夠大的浮力可以引起裂紋向上傳播.Secor和Pollard（1975）研究了浮力驅(qū)動斷裂作為巖漿傳播的機制.Spence和Sharp（1985）得到了一個無窮大板中流體源作用下流體驅(qū)動斷裂的相似解，但沒有考慮浮力.Spence等（1987）推導出一個彈性固體中半無窮含流體裂紋浮力驅(qū)動斷裂的穩(wěn)態(tài)解，他們使用了潤滑理論描述流體在斷裂中的黏性流動，并應(yīng)用微分方程描述彈性變形.Lister（1990）考慮流體恒常流量下浮力驅(qū)動裂紋傳播的材料韌性效應(yīng)，得到了一個裂紋表面位移解，并表明巖脈為泡狀頭和等寬度長尾巴的形狀.Lister和Kerr（1991）描述了巖脈傳播中流體和固體耦合效應(yīng)，指出巖脈中巖漿傳輸是受浮力和黏性壓力降之間的平衡所控制，除了在巖脈尖端外，彈性力起著次級作用.Rubin（1993）得到了一個在給定壓力源下巖脈增長的自相似解.Rubin（1995）綜述了對圍巖破裂包括觀察、物理過程和未來方向等各個方面的內(nèi)容.Meriaux和Jaupart（1998）考慮有限彈性板中浮力和給定超壓力驅(qū)動下巖脈傳播.Bonafede和Rivalta（1999）使用位錯模型研究巖脈周圍的應(yīng)力和位移場.Roper和Lister（2005）進一步研究含流體裂紋在超高壓半無窮彈性板的傳播過程，得到了一個臨界長度，即當巖脈長度小于這個臨界長度時，巖漿房超高壓起著主要作用，當巖脈長度大于這個臨界長度時，浮力起著主要作用.Chen等（2007）利用攝動法計算了地殼密度逐漸變化條件下巖脈穩(wěn)態(tài)傳播過程，并分析了裂紋穩(wěn)態(tài)傳播的條件.

前人有利用彈性半空間的點源模型來模擬研究火山區(qū)巖漿囊膨脹與地表變形之間的關(guān)系（Mogi，1958），有利用黏彈性圍巖球殼模型研究球狀巖漿囊周圍應(yīng)力應(yīng)變響應(yīng)（Dragoni and Magnanensi，1989），及利用黏彈性扁橢球模型研究巖漿囊膨脹對斷層的影響（Newman etal.，2006）.本文所研究的巖漿沿已有斷層運移主要涉及黏性巖漿流體在壓應(yīng)力條件下打開已有斷層并在彈性地殼中上升的過程，我們采用三維彈性上地殼模型進行模擬.

已有分析表明，在構(gòu)造應(yīng)力和巖石的自重作用下，殼內(nèi)以分異演化后的中酸性巖漿為主的巖漿與圍巖之間存在密度差，巖漿會受到一個向上的浮力.由此可能導致上覆圍巖中已有斷層受到接觸巖漿的擠壓應(yīng)力時，斷層張開，巖漿向上運移.在此過程中，靜巖壓力、構(gòu)造應(yīng)力以及密度差造成的浮力會隨之減小，有可能造成巖漿運移的終止.如果巖漿受到的超壓力足夠大，則巖漿能夠向上運移直至地表.

如上所述，過去關(guān)于巖漿在地殼中運移的研究多采用斷裂力學方法，研究巖漿驅(qū)動下圍巖破裂的巖漿傳播過程，且大多是二維的.本文針對巖漿在已有斷層中運移過程，采用三維有限元法與有限差分法相結(jié)合，計算模擬巖漿在地殼中的傳播過程.我們利用自行編制的相關(guān)程序，對長白山天池火山進行了參數(shù)試驗的模擬計算.有限差分方法用來分析斷層中巖漿的流動，有限元方法用來計算巖脈傳播過程中已有斷層的打開導致巖體變形過程.

2 物理模型與研究方法

巖漿沿著已知斷層面向上運移便形成了巖脈.本文用兩個相鄰的正方體模型來模擬巖體，相鄰的接觸面作為斷層面，模擬巖漿沿已有斷層運移的過程.在水平方向上巖體受到周圍構(gòu)造應(yīng)力的作用；巖漿在斷層中受到構(gòu)造應(yīng)力與彈性形變產(chǎn)生的壓應(yīng)力，以及靜巖壓力與自身重力.在三維模型下，巖漿在斷層中運移的寬度（沿斷層面）是有限且不變的.

考慮從巖漿房中垂直向上進入圍巖的巖漿，假設(shè)巖漿為不可壓縮的黏性流體，且在運移過程中不會發(fā)生固化，利用流體力學的潤滑理論描述侵入斷層中的巖漿流體（Rice，1968），在流體充滿的空間中，根據(jù)質(zhì)量守恒，有

這里w是斷層面張開位移，q是流體流量.上式寫成差分格式為

其中wi-w0i是第i個斷層單元當前張開位移與Δt時間前張開位移之差，Δqi是第i個單元的流量差分.

根據(jù)流體流量與孔隙壓力的Poiseuille′s定律，有

其中，μ是巖漿的動力黏度，q是巖脈中巖漿流量，w是巖漿侵入使斷層張開的位移，p是巖脈中巖漿受到的壓應(yīng)力，ρm是巖漿密度，z表示沿巖脈向上運移方向.驅(qū)動巖脈傳播壓力，即為使巖石彈性形變產(chǎn)生的壓應(yīng)力，包括巖漿流體的壓應(yīng)力、巖石圍壓和構(gòu)造應(yīng)力，因此式（3）中

其中pe＝pe（z）為巖漿壓力，使圍巖產(chǎn)生彈性形變，ρs是圍巖密度，Δσ為底部垂直于巖脈平面的構(gòu)造應(yīng)力（Rubin，1995；Roper and Lister，2005）.（4）式可寫為

其中

（3）式成為

其中

由（7）式，有

故

同樣，

其中pe的導數(shù)寫為差分格式為

由（9）—（12）得q 的中心差分格式為

結(jié)合（2）式，得

其中

對于（14）式，考慮彈性壓力pe（z）與張開位移w（z）之間的彈性關(guān)系，當給定一組張開位移值時，可以確定一組彈性壓力.其中的w初始給定為0.005m，pe兩端的條件為

其中S0為侵入巖漿底部構(gòu)造應(yīng)力與靜巖壓力之和，ΔP為給定的超壓力，來自于巖漿流體的內(nèi)部壓力與浮力，本研究采用恒常的超壓值.S（m）為隨運移過程變化的侵入巖漿尖端構(gòu)造應(yīng)力與靜巖壓之和.

對于非連續(xù)變形分析中的接觸問題，當pe（z）已知，w（z）可根據(jù)有限單元方法（王勖成，2003）求出.假設(shè)地殼中存在一斷層，當侵入巖漿末端的壓力pe（a）大于靜巖壓與構(gòu)造應(yīng)力之和時，已存在的斷層會張開，巖漿便會繼續(xù)侵入并向上運移.具體計算時，先由已有斷層初始的微小張開位移w0，利用有限差分格式計算出平衡公式（14）中的應(yīng)力項，得到巖漿壓力pe，從而獲得圍巖的應(yīng)力狀態(tài).得到圍巖-巖漿接觸面上的邊界應(yīng)力后，用有限單元方法得到此應(yīng)力邊界條件下圍巖中各單元的應(yīng)變量，進而獲得接觸面上的位移值，即已有斷層兩側(cè)的張開位移量wi.如果每次獲得的接觸邊界上最上端第m個單元端部節(jié)點的張開位移值wm+j大于零，則可以判斷斷層繼續(xù)向上張開，否則，巖漿運移停止.如果斷層進一步張開，巖漿侵入，則以新獲得的wi作為下次計算中的w0，重復上述計算，直到斷層面停止張開為止.采用有限單元法計算接觸面的位移值，相當于計算巖脈前端巖石所受應(yīng)力.當其受拉應(yīng)力作用時斷層張開，巖脈向上運移.上述過程將包含已有斷層的彈性問題轉(zhuǎn)化為了一個接觸問題.

3 計算模型

考慮一個含有已知斷層的區(qū)域巖體，構(gòu)建一個三維長方體區(qū)域，底面為共用一條邊的兩個正方形，正方形邊長為10km.長方體高h（h＝5000、8000、10000m），構(gòu)成相互接觸的兩個正方體區(qū)域（圖1a）.如圖劃分網(wǎng)格單元，x方向（斷層面方向）中心適當加密，y方向均勻分割，z方向固定為49等份（50層節(jié)點），共計3528個單元和5000個單元節(jié)點.兩正方體的接觸面模擬為斷層面，在區(qū)域下方有一巖漿房，設(shè)置初始進入斷層微小張開的巖漿位于長方體的底面中心單元并向上延伸5層節(jié)點單元.三維模型下巖漿充滿斷層一個單元的寬度（圖1中沿x軸方向）.

圖1 （a）計算模型網(wǎng)格，底部紅色代表侵入斷層的巖漿；（b）侵入巖漿三維示意圖，w為巖漿侵入使斷層兩側(cè)張開的位移；（c）模型受構(gòu)造應(yīng)力（Δσ）示意圖，p0為底部節(jié)點的巖漿壓力；（d）侵入巖漿二維示意圖，w為斷層兩側(cè)張開位移，pe、p0為巖漿壓力，a為巖漿侵入高度Fig.1 （a）Grid model for simulation with red bar for invaded magma；（b）Three dimension schematic for invaded magma，wstand for open displacement of the crack；（c）Model schematic with tectonic stress（Δσ）；（d）Two dimension schematic for invaded magma，wstand for open displacement of the crack，peand p0for magma pressure，afor invasion height of magma

模型的底面節(jié)點垂直方向固定，四周各面加載10MPa左右的構(gòu)造應(yīng)力，且x方向上略大于y方向.構(gòu)造應(yīng)力與靜巖壓力之合力沿垂直方向線性變化.斷層從下部開始張開5層節(jié)點的初始微小位移w0＝0.005m，最下層處節(jié)點的巖漿壓力為p0，見式（16）.超壓力ΔP為固定大小，當模型深度h一定時，p0恒定不變.

本文所采用的圍巖彈性模量E＝50GPa，泊松比ν＝0.25，重力加速度g＝9.8m·s-2.我們對不同組參數(shù)（包括巖漿房超壓力，巖漿與圍巖的密度差及構(gòu)造應(yīng)力）組合進行了數(shù)值計算.

4 計算結(jié)果與討論

在原有三維模型的基礎(chǔ)上，改變設(shè)置得到二維模型并進行相關(guān)計算，用以對比三維的情形.二維模型采用類似平面應(yīng)變的模型，僅僅將三維模型中一個單元內(nèi)的巖漿寬度擴展為沿x方向（斷層面走向）上充滿斷層面，其他條件相似.

兩個算例模型h＝8km，σx＝80MPa，σy＝73.2MPa，Δρ＝500kg·m-3，μ＝2.0Pa·s，ΔP＝17MPa.

3-D模型中巖脈的張開位移為3～5m，向上運移逐漸呈現(xiàn)火把形截面；2-D模型張開位移為10～120m，截面形狀為盾形或梯形（圖2）.兩個模型中，巖脈在接近地表時上升速度都會加快，這是因為淺層靜巖壓力和構(gòu)造應(yīng)力均減小之故.2-D模型計算得到的張開位移量與根據(jù)斷裂力學二維計算的結(jié)果一致（Chen etal.，2007）.

表面（地表）垂直位移的形態(tài)上，3-D模型先由斷層兩側(cè)緩慢隆起，進而中心漸漸凹陷，兩側(cè)繼續(xù)突出，使斷層與斷層兩側(cè)突出之間高差增大，最終斷層中心較快速突出，形成山峰狀；2-D模型則先由斷層兩側(cè)緩慢抬升，進而中間凹陷成低地，兩側(cè)繼續(xù)抬升，使凹陷與兩側(cè)抬升之間高差增大，最終在斷層上方較快速地突出，形成山脊狀.三維模型中地表最大垂直位移值在0.3m以內(nèi)；二維模型中地表最大垂直位移值達到12m（脊）和20m（槽）（圖3）.二維模型中的地表位移總體上明顯大于三維模型中的，且在前期位移量均較?。?-D中為mm～cm級，2-D中為dm～m級），在巖脈接近地表時垂直位移量增大.二維模型中的垂直位移比三維模型的垂直位移大約兩個數(shù)量級.

兩種模型的pe值曲線基本相同（圖4），均是先穩(wěn)定增大，在接近巖脈頂端附近迅速減小至S（m），即構(gòu)造應(yīng)力與靜巖壓力之和（見式（16））.這是由于在每一步巖脈運移之后給定的應(yīng)力條件限制的.

三維模型中如果限定侵入斷層的巖漿寬度適當小時，可以模擬斷裂帶的火山巖漿運移過程.二維模型可以近似模擬洋中脊巖漿上涌的過程.

我們進行了不同參數(shù)下的三維模擬計算，研究了巖漿黏度、密度差、模型深度對最小超壓力（巖漿運移到地表所需的最小巖漿超壓力）的影響.結(jié)果表明，對于深度為10km的模型，相應(yīng)底部邊界的圍壓（構(gòu)造應(yīng)力與靜巖壓力）為σx＝100MPa，σy＝91.5MPa，巖漿的動力黏度取為2×10-7～2×10-3MPa·s.當超壓力達到17～20MPa時，可以驅(qū)動巖脈運移到地表.同時，巖漿動力黏度越大，使巖脈運移到地表需要的最小超壓力就越大（見表1）.隨著模型深度的增加，黏度增大對最小超壓力的影響減弱，即增大黏度后，所需最小超壓力隨模型深度增加緩慢，且在一定黏度范圍內(nèi)，最小超壓力變化不大.密度差（圍巖有效密度與巖漿密度之差，見式（8））取值為300～700kg·m-3，影響趨勢為密度差越小，所需最小超壓力越大，且影響幅度比較小——小于1MPa／（0.1g·cm-3）.

對于深度小于10km的模型進行的計算表明，當巖漿囊埋藏較淺時，使巖脈發(fā)生運移并運移到地表的最小超壓力要比上述10km情形下有明顯減小.也許是因為較淺部構(gòu)造應(yīng)力和靜巖壓力的影響減小，使得開始發(fā)生運移及維持運移都較為容易.

5 長白山天池火山的應(yīng)用

長白山天池火山是長白山火山區(qū)最具代表性的三座巨型堿性火山之一，是該區(qū)發(fā)育最為完全的火山之一，是我國境內(nèi)保存最好的新生代多成因中央式火山（魏海泉，2010）.天池火山從中晚更新世至今有過多次噴發(fā)，最近一次大噴發(fā)在1215年前后，其規(guī)模和程度被認為是近2000年來地球上最大的爆炸式火山噴發(fā)之一.自1985年開展天池火山的地震監(jiān)測以來，2002年7月以前長白山天池火山地區(qū)火山地震活動性一直比較低，從2002年7月開始到2004年9月，該地區(qū)地震活動性明顯增強，地形變也出現(xiàn)了較為明顯的變化，這些現(xiàn)象都表明長白山天池火山地區(qū)地下巖漿或熱液活動正在發(fā)生變化（明躍紅等，2006），相應(yīng)于天池火山地下淺層殼內(nèi)巖漿囊的壓力與體積變化.

表1 不同參數(shù)下使巖漿運移到地表的最小超壓力Table 1 The least overpressure to drive magma up to surface under different parameters

前人對長白山火山區(qū)的地質(zhì)結(jié)構(gòu)進行了深入的研究，對于該地區(qū)地殼內(nèi)巖漿囊的賦存有了進一步認識（劉若新等，1998a；樊祺誠等，2005；張成科等，2002；湯吉等，2001；楊卓欣等，2005；趙金仁等，2003）.近年來的一些研究顯示出該地區(qū)巖漿囊的新特點.張先康等（2002）利用人工地震測深得到的分層巖漿系統(tǒng)最淺位于8～9km；吳建平等（2007）對天池火山區(qū)震群活動的研究指出，5km深度附近可能存在巖漿熱液活動或與巖漿增壓有關(guān)；朱桂芝等（2008）利用Mogi模型聯(lián)合GPS觀測數(shù)據(jù)并結(jié)合遺傳算法反演得到長白山火山區(qū)巖漿囊位于深度9.2 km處；陳國滸等（2008）基于InSAR、GPS形變場的Mogi模型反演得到的結(jié)果推測出天池火山下方7.9km存在一個巖漿壓力源；吳建平等（2009）利用面波層析成像得到天池火山附近S波低速層頂層埋深在8～9km；魏海泉（2010）參考有花崗質(zhì)巖石提出的模型，指出殘余的粗面質(zhì)熔體在5km或更深處的地殼內(nèi)聚集.這些研究顯示出長白山火山區(qū)地殼淺部的巖漿囊或許存在于5～10km處.故此我們分別以h＝10km和h＝5km的模型來模擬天池火山的巖漿運移并進行相關(guān)討論.

天池火山巖漿是經(jīng)過部分礦物結(jié)晶分異演化的巖漿，因此其幔源的性質(zhì)減弱，而偏向于中、酸性的粗面巖質(zhì)及堿流巖質(zhì)類型（劉若新等，1998b）.從火山活動的噴出物來看，由早期（上新世到早更新世）的造盾，到后期（早更新世到晚更新世）的造錐，再到全新世-近代以來的浮巖噴發(fā)，逐漸由偏基性的玄武質(zhì)巖石為主過渡到偏中、酸性（粗面巖、堿流巖）的粗面-堿流質(zhì)巖石為主（劉若新等，1998a；劉若新等，1998b；李春鋒等，2006）.從巖石化學組成上看，早期造盾時期以高鎂鐵性質(zhì)的拉斑系列巖石為主；后期造錐時期，以堿性、鈣堿性的粗面巖、流紋巖為主.浮巖的噴發(fā)也以粗面質(zhì)和堿流質(zhì)居多（劉若新等，1998a；樊祺誠等，2005）.巖漿的分異演化伴隨著各個時期的火山活動.天池火山巖漿也存在著混合作用.來自年代學的研究顯示，粗面玄武巖的噴發(fā)從早期造盾階段結(jié)束后，一直持續(xù)到近代的噴發(fā).對來自該區(qū)域的巖石及深淺條帶混合等現(xiàn)象的分析表明（樊祺誠等，2005），天池火山造錐噴發(fā)的不同時期，存在著粗面質(zhì)巖漿混合現(xiàn)象.有人認為千年大噴發(fā)的觸發(fā)機制可能是深部的地幔玄武巖漿注入地殼巖漿房的結(jié)果（劉若新等，1998b）.水含量的多少對巖漿黏度有顯著影響.劉若新等研究發(fā)現(xiàn)千年大噴發(fā)中堿流質(zhì)巖漿中的水含量已達到飽和甚至過飽和（劉若新等，1998b），這些因素都會使淺部巖漿的黏度變得更?。‵reundt and Schmincke，1992；Hui，2008）.考慮粗面質(zhì)巖漿的混合以及可能的不同水含量影響，巖漿黏度值可達到104～109Pa·s.

圖2 巖漿運移上升圖左側(cè)：3-D模型結(jié)果，右側(cè)：2-D模型結(jié)果；由下至上：巖漿運移隨時間變化截圖，橫軸為y向（見圖1a），縱軸為z向.Fig.2 Magma migration up to surface Left：3-D model result，right：2-D model result；Bottom to up：time screenshots of migration with horizontal yand vertical z（see Fig.1a）.

在用本文上述模型模擬長白山天池火山的運移時，取2×106Pa·s作為模型中巖漿的黏度值.本文就巖漿房可能在深度10km和5km兩種情形計算了巖漿發(fā)生運移達到地表時所需要的最小超壓力，同時模擬了巖漿運移過程中地表隨時間的變化特征，并根據(jù)巖漿房壓力可能由于巖漿流出而導致減小的理論關(guān)系，估算了相應(yīng)情形下使巖漿運移至地表時地下巖漿囊所需的體積大小.

圖3 地表垂直位移變化圖左側(cè)：3-D模型結(jié)果，右側(cè)：2-D模型結(jié)果；由下至上：地表垂直位移隨時間變化截圖，色標為位移值.Fig.3 Vertical displacement variation on surface during magma migration Left：3-D model result，right：2-D model result；Bottom to up：time screenshots of migration with vertical displacement value u.

若使巖漿運移到地表，巖漿囊頂部的超壓力需要達到21MPa（深度10km）和10MPa（深度5km），可以看出較淺模型需要的最小超壓力要小.前人得到的天池火山千年大噴發(fā)時巖漿房超壓為6.25MPa和15MPa（魏海泉，2010），前者對應(yīng)于5km深的堿流質(zhì)層狀巖漿房，后者對應(yīng)粗面質(zhì)層狀巖漿房（應(yīng)位于更深處）.較淺部巖漿囊受到的構(gòu)造應(yīng)力較小，達到運移條件所需的最小超壓力也較小.

地表位移變化與上述模型中給出的結(jié)果相似（參考圖3-左圖）.由地表位移變化圖可以看出，地表的隆起隨著巖漿向上運移而發(fā)生變化，仍然表現(xiàn)出斷層面對應(yīng)巖脈兩側(cè)先緩慢隆起，而后向巖脈中心聚攏.地表垂直位移沒有呈現(xiàn)中心放射狀，原因首先在于考慮到模型的簡化以及計算的條件，本文模型設(shè)置單一已有斷層，而長白山天池火山處于該地區(qū)兩大主要斷裂——北東向的六道溝—天池—甑峰山斷裂和北西向的白山鎮(zhèn)—天池—金策斷裂——的交叉點，這兩條走向近乎垂直的斷裂控制著天池火山的產(chǎn)出（劉明軍等，2004）；其次模型上部沒有預設(shè)破火山口狀的隆起，而對于天池火山已有中心隆起的地形，這類效應(yīng)或許不會很大.從巖漿運移到地表后的最終形態(tài)來看，發(fā)生較大垂直位移的范圍大致呈現(xiàn)橢圓的形態(tài)，也與上述原因有關(guān).

圖4 巖漿壓力pe變化圖左側(cè)：3-D模型結(jié)果，右側(cè)：2-D模型結(jié)果；由下至上：pe隨巖漿上升變化，z為巖漿上升高度.Fig.4 The variation of magma pressure pe Left：3-D model result，right：2-D model result.Bottom to up：time screenshots of pevariation during magma migration with magma height z.

在本文模型的條件下，排除大規(guī)模構(gòu)造運動使得構(gòu)造應(yīng)力發(fā)生劇烈變化的情況，若給定適當?shù)某瑝毫?，則可以粗略估算這種運移發(fā)生時天池火山下方淺部巖漿囊的大小.可以利用下列公式（Maccaferri etal.，2010）估算模型下方巖漿囊的容量變化與超壓力變化的關(guān)系，

其中V為巖漿囊中的巖漿體積V0與侵入巖脈中的巖漿體積Vout之和，Kp為巖漿的體積模量，ΔP為超壓力.上式可改寫為

根據(jù)上式將模型中的超壓力設(shè)置為隨著巖漿向上運移而發(fā)生變化.對于不同深度的模型取最小超壓力作為初始ΔP0，即5km時ΔP0＝10MPa；10km時ΔP0＝21MPa.巖漿體積模量Kp＝10000MPa.

由每一步巖脈向上運移的高度可算出該次運移由巖漿房向已有斷層中侵入的巖漿體積，即為Vout.為保證巖漿可以在給定初始超壓力下能夠運移至地表，計算出應(yīng)給定的在運移發(fā)生前巖漿囊的體積.對于5km的模型，如果巖漿囊超壓力為10MPa，下方巖漿囊體積大于6×109m3時可以保證巖漿運移到地表；對于10km的模型，如果巖漿囊超壓力為21MPa，下方巖漿囊容量為2×1010m3時可以保證巖漿運移到地表.如果考慮層狀巖漿囊，假設(shè)半徑為2～6km，則此時5km模型巖漿囊的厚度約為50～500m，10km模型巖漿囊厚度約為200～1600m.

上述對巖漿囊大小的估計，是在構(gòu)造應(yīng)力不發(fā)生劇烈變化且給定超壓力情況下所作的計算，如果巖漿能夠運移到地表，那么天池火山下部對應(yīng)深度可能存在的巖漿囊原有的大小.

6 結(jié)論

對于10km深度的巖漿囊，考慮巖漿浮力和構(gòu)造應(yīng)力情況下，巖漿黏度在2×10-3到2×10-7MPa·s范圍時（相當于玄武巖漿），當超壓力達到17～20MPa時，巖漿可以沿著已有斷層向上運移至地表，從而形成巖脈.在總體上看，巖漿黏度越大，巖脈發(fā)生運移且至地表的最小超壓力就越大；密度差（圍巖有效密度和巖漿密度之差）越小，所需最小超壓力就越大.

巖漿囊埋藏越淺（10～8km），使巖脈發(fā)生運移且至地表的最小超壓力就越小.

三維模型中巖脈的張開位移為3～5m，二維模型中則為10～120m.在接近地表時，兩種模型中巖脈的運移速度都會加快.

三維模型得到的地表垂直位移，由中心凹陷、兩側(cè)突出逐漸形成中心突出的峰狀；二維模型中則由中間凹陷為槽狀逐漸形成中間凸起的脊狀.三維模型可以近似模擬斷裂帶控制的火山巖漿活動，二維模型可以近似模擬洋中脊的情況.

用本文模型模擬長白山天池火山巖漿運移時，地表垂直位移沒有呈現(xiàn)中心放射狀，原因在于設(shè)置的單一已有斷層，而天池火山處于該地區(qū)兩大斷裂的交叉點，這兩條走向近乎垂直的斷裂控制著天池火山的產(chǎn)出；其次模型上部沒有預設(shè)破火山口狀的隆起，而對于天池火山已有中心隆起的地形，這類效應(yīng)或許不會很大.

若假設(shè)長白山天池火山下5km或者10km深度存在巖漿囊，應(yīng)用本文的模型，如果巖漿要運移至地表，超壓力至少應(yīng)為10MPa和21MPa.當構(gòu)造應(yīng)力不發(fā)生劇烈變化且給定超壓力分別為10MPa和21MPa時，若巖漿運移到地表，巖漿囊應(yīng)有的最小體積約為6×109m3（5km深）和2×1010m3（10km深）.相對應(yīng)于半徑為2～6km的層狀巖漿囊時，厚度約為50～500m（5km深）和200～1600m（10km深）.

本文對于長白山天池火山地區(qū)巖漿活動機制和特征及噴發(fā)條件的研究結(jié)果為未來該地區(qū)火山噴發(fā)預測和監(jiān)控提供了參考.

致謝 感謝中國科學院地球深部研究重點實驗室提供計算模擬平臺.感謝中國科學院超級計算中心提供科學計算網(wǎng)格服務(wù).

Bonafede M，Rivalta E.1999.On tensile cracks close to and across the interface between two welded elastic half-space.Geophys.J.Int.，138（2）：410-434，doi：10.1046／j.1365-246X.1999.00880.x.

Chen G H，Shan X J，Moon W M，etal.2008.A modeling of the magma chamber beneath the Changbai Mountains volcanic area constrained by InSAR and GPS derived deformation.Chinese J.Geophys.（in Chinese），51（4）：1085-1092，doi：10.3321／j.issn：0001-5733.2008.04.017.

Chen Z，Jin Z H，Johnson S E.2007.A perturbation solution for dyke propagation in an elastic medium with graded density.Geophys.J.Int.，169（1）：348-356，doi：10.1111／j.1365-246X.2006.03297.x.

Dragoni M，Magnanensi C.1989.Displacement and stress produced by apressurized，Spherical magma chamber，Surrounded by a viscoelastic shell.Physics of the Earth and Planetary Interiors，56（3-4）：316-328，doi：10.1016／0031-9201（89）90166-0.

Fan Q C，Sui J L，Sun Q，etal.2005.Preliminary research of magma mixing and explosive mechanism of the Millennium eruption of Tianchi volcano.Acta Petrologica Sinica（in Chinese），21（6）：1703-1708，doi：10.3321／j.issn：1000-0569.2005.06.017.

Freundt A，Schmincke H U.1992.Mixing of rhyolite，trachyte and basalt magma erupted from a vertically and laterally zoned reservoir，composit flow P1，Gran Canaria.Contrib.Mineral.Petrol.，112（1）：1-19，doi：10.1007／BF00310952.

Hui H J.2008.Viscosity of Silicate Melts［Ph.D.thesis］.Ann Arbor，USA：The University of Michigan.

Li C F，Zhang X K，Zhang Y，etal.2006.Analysis of tectonic setting of Changbaishan Tianchi volcano.Seismological and Geomagnetic Observation and Research（in Chinese），27（5）：43-49，doi：10.3969／j.issn.1003-3246.2006.05.009.

Lister J R.1990.Buoyancy-driven fluid fracture：the effects of material toughness and of low-viscosity precursors.J.Fluid Mech.，210：263-280，doi：10.1017／S0022112090001288.

Lister J R，Kerr R C.1991.Fluid-mechanical models of crack propagation and their application to magma transport in dykes.J.Geophys.Res.，96（B6）：10049-10077，doi：10.1029／91JB00600.

Liu M J，Gu M L，Sun Z G，etal.2004.Activity of main faults and hydrothermal alteration zone at the Tianchi volcano，Changbaishan.Earthquake Research in China（in Chinese），20（1）：64-72，doi：10.3969／j.issn.1001-4683.2004.01.007.

Liu R X，F(xiàn)an Q C，Zheng X S，etal.1998a.Magma evolution of Changbaishan volcano.Science in China（Series D） （in Chinese），28（3）：226-231

Liu R X，Wei H Q，Li J T，etal.1998b.The Recent Eruption of Tianchi Volcano，Changbiashan（in Chinese）.Beijing：Science Press.

Maccaferri F，Bonafede M，Rivalta E.2010.A numerical model of dyke propagation in layered elastic media.Geophys.J.Int.，180（3）：1107-1123，doi：10.1111／j.1365-246X.2009.04495.x.

Meriaux C，Jaupart C.1998.Dike propagation through an elastic plate.J.Geophys.Res.，103（B8）：18295-18314，doi：10.1029／98JB00905.

Ming Y H，Su W，F(xiàn)ang L H.2006.A preliminary study of types of volcanic earthquakes and volcanic activity in Changbaishan，Tianchi volcano.Earthquake Research in China（in Chinese），22（1）：56-63，doi：10.3969／j.issn.1001-4683.2006.01.006.

Mogi K.1958.Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them.Bulletin of the Earthquake Research Institute，36：99-134.

Newman A V，Dixon T H，Gourmelen N.2006.A four-dimensional viscoelastic deformation model for long valley caldera，california，between 1995and 2000.Journal of Volcanology and Geothermal Research，150（1-3）：244-269，doi：10.1016／j.jvolgeores.2005.07.017.

Rice J R.1968.Mathematics analysis in the mechanics of fracture.／／Liebowitz H.Fracture，an Advanced Treatise，Vol.Ⅱ，Chapter 3.New York：Academic Press，191-311.

Roper S M，Lister J R.2005.Buoyancy-driven crack propagation from an overpressure source.J.Fluid Mech.，536：79-98，doi：10.1017／S0022112005004337.

Rubin A M.1993.Dykes vs.diapirs in viscoelastic rock.Earth Planet.Sci.Lett.，119（4）：641-659，doi：10.1016／0012-821X（93）90109-M.

Rubin A M.1995.Propagation of magma-filled crack.Annu.Rev.Earth Planet.Sci.，23：287-333，doi：10.1146／annurev.ea.23.050195.001443.

Secor D，Pollard D.1975.On the stability of open hydraulic fractures in the earth′s crust.Geophys.Res.Lett.，2：510-513，doi：10.1029／GL002i011p00510.

Spence D A，Sharp P.1985.Self-similar solutions for elastohydrodynamic cavity flow.Proc.R.Soc.Lond.，400（189）：289-313，doi：10.1098／rspa.1985.0081.

Spence D A，Sharp P W，Turcotte D L.1987.Buoyancy-driven crack propagation：a mechanism for magma migration.J.Fluid Mech.，174：135-153，doi：10.1017／S0022112087000077.

Tang J，Deng Q H，Zhao G Z，etal.2001.Electric conductivity and magma chamber at the Tianchi volcano area in Changbaishan mountain.Seismology and Geology（in Chinese），23（2）：191-200，doi：10.3969／j.issn.0253-4967.2001.02.008.

Wang X C.2003.Finite Element Method（in Chinese）.Beijing：Tsinghua University Press.

Weertman J.1971.Theory of water-filled crevasses in glaciers applied to vertical magma transport beneath oceanic ridges.J.Geophys.Res.，76（5）：1171-1183，doi：10.1029／JB076i005p01171.

Wei H Q.2010.Magma up-moving process within the magma prism beneath the Changbaishan volcanoes.Earth Science Frontiers（in Chinese），17（1）：011-023.

Wu J P，Ming Y H，Zhang H R，etal.2007.Earthquake swarm activity in Changbaishan Tianchi volcano.Chinese J.Geophys.（in Chinese），50（4）：1089-1096，doi：10.3321／j.issn：0001-5733.2007.04.016.

Wu J P，Ming Y H，Su W，etal.2009.S wave velocity structure beneath Changbaishan Tianchi volcano inferred from receiver function.Seismology and Geology（in Chinese），31（4）：584-597，doi：10.3969／j.issn.0253-4967.2009.04.002.

Yang Z X，Zhang X K，Zhao J R，etal.2005.Tomographic imaging of 3-D crustal structure beneath Changbaishan-Tianchi Region.Chinese J.Geophys.（in Chinese），48（1）：107-115，doi：10.3321／j.issn：0001-5733.2005.01.016.

Zhang C K，Zhang X K，Zhao J R，etal.2002.Study on the crustal and upper mantle structure in the Tianchi volcanic region and its adjacent area of Changbaishan.Chinese J.Geophys.（in Chinese），45（6）：812-820，doi：10.3321／j.issn：0001-5733.2002.06.008.

Zhang X K，Zhang C K，Zhao J R，etal.2002.Deep seismic sounding investigation into the deep structure of the magma system in Changbaishan-Tianchi volcano region.Acta Seismologica Sinica（in Chinese），15（2）：135-143，doi：10.3321／j.issn：0253-3782.2002.02.003.

Zhao J R.2003.3Dtomography of velocity structure in upper crust beneath the Changbaishan Tianchi volcanic region.Chinese J.Geophys.（in Chinese），46（6）：796-802，doi：10.3321／j.issn：0001-5733.2003.06.011.

Zhu G Z，Wang Q L，Shi Y L，etal.2008.Modelling pressurized deformation source for Changbaishan volcano with homogenous expansion point source.Chinese J.Geophys.（in Chinese），51（1）：108-115，doi：10.3321／j.issn：0001-5733.2008.01.014.

附中文參考文獻

陳國滸，單新建，Moon W M等.2008.基于InSAR、GPS形變場的長白山地區(qū)火山巖漿囊參數(shù)模擬研究.地球物理學報，51（4）：1085-1092，doi：10.3321／j.issn：0001-5733.2008.04.017.

樊祺誠，隋建立，孫謙等.2005.天池火山千年大噴發(fā)的巖漿混合作用與噴發(fā)機制初步探討.巖石學報，21（6）：1703-1708，doi：10.3321／j.issn：1000-0569.2005.06.017.

李春鋒，張興科，張旸等.2006.長白山天池火山的地質(zhì)構(gòu)造背景.地震地磁觀測與研究，27（5）：43-49，doi：10.3969／j.issn.1003-3246.2006.05.009.

劉明軍，顧夢林，孫振國等.2004.長白山天池火山主斷裂活動與熱液蝕變帶.中國地震，20（1）：64-72，doi：10.3969／j.issn.1001-4683.2004.01.007.

劉若新，樊祺誠，鄭祥身等.1998a.長白山天池火山的巖漿演化.中國科學（D輯），28（3）：226-231.

劉若新，魏海泉，李繼泰等.1998b.長白山天池火山近代噴發(fā).北京：科學出版社.

明躍紅，蘇偉，房立華.2006.長白山天池火山地震類型及火山活動性的初步研究.中國地震，22（1）：56～63，doi：10.3969／j.issn.1001-4683.2006.01.006.

湯吉，鄧前輝，趙國澤等.2001.長白山天池火山區(qū)電性結(jié)構(gòu)和巖漿系統(tǒng).地震地質(zhì)，23（2）：191-200，doi：10.3969／j.issn.0253-4967.2001.02.008.

王勖成.2003.有限單元法.北京：清華大學出版社.

魏海泉.2010.長白山火山巖漿柱巖漿上升作用過程.地學前緣，17（1）：11-23.

吳建平，明躍紅，張恒榮等.2007.長白山天池火山區(qū)的震群活動研究.地球物理學報，50（4）：1089-1096，doi：10.3321／j.issn：0001-5733.2007.04.016.

吳建平，明躍紅，蘇偉等.2009.長白山火山區(qū)殼幔S波速度結(jié)構(gòu)研究.地震地質(zhì)，31（4）：584-597，doi：10.3969／j.issn.0253-4967.2009.04.002.

楊卓欣，張先康，趙金仁等.2005.長白山天池火山區(qū)三維地殼結(jié)構(gòu)層析成像.地球物理學報，48（1）：107-115，doi：10.3321／j.issn：0001-5733.2005.01.016.

張成科，張先康，趙金仁等.2002.長白山天池火山區(qū)及鄰近地區(qū)殼幔結(jié)構(gòu)探測研究.地球物理學報，45（6）：812-820，doi：10.3321／j.issn：0001-5733.2002.06.008.

張先康，張成科，趙金仁等.2002.長白山天池火山區(qū)巖漿系統(tǒng)深部結(jié)構(gòu)的深地震測深研究.地震學報，15（2）：135-143，doi：10.3321／j.issn：0253-3782.2002.02.003.

趙金仁.2003.長白山天池火山區(qū)上地殼三維速度層析成像.地球物理學報，46（6）：796-802，doi：10.3321／j.issn：0001-5733.2003.06.011.

朱桂芝，王慶良，石耀霖等.2008.各向同性膨脹點源模擬長白山火山區(qū)巖漿囊壓力變形源.地球物理學報，51（1）：108-115，doi：10.3321／j.issn：0001-5733.2008.01.014.