地壳中岩浆沿已有断层运移的数值模拟——在长白山天池火山的应用

程旭，陈祖安，白武明

1中国科学院地质与地球物理研究所地球深部研究重点实验室，北京 100029

2中国科学院大学，北京 100049

1 引言

幔源的岩浆在上升侵入地壳之后，会在不同深度形成岩浆囊或岩浆房.这些岩浆囊与下部的深层岩浆房仍有通道连接.当不断有深部岩浆补给到浅部的岩浆囊中时，岩浆囊中的压力就会持续增大，由此可能会导致岩浆运移.岩浆在地壳中的运移分为两种类型：一种是岩浆压力驱动下围岩不断破裂，岩脉不断增长；另一种是岩浆沿着已有断层或通道运移.

岩浆运移机制及与之相关的围岩破裂是一个岩浆从开始缓慢侵入逐渐过渡到急速上涌喷发的过程.要全面认识这一过程，除了需要先进的仪器来监测深部应力外，解析分析与数值模拟流体与固体破裂相耦合的这一动力学过程也是相当重要的.关于岩浆在地壳中运移过程的力学已经取得了许多研究成果.关于从岩浆囊开始的岩脉传播的研究大部分集中在线弹性断裂力学框架上的断裂扩展过程，这些同时也考虑了岩浆流动的流体动力学过程.Weertman（1971）考虑了一个受到张应力和重力诱发的静水压力下的水平弹性板中含流体裂纹的问题，研究结果表明足够大的浮力可以引起裂纹向上传播.Secor和Pollard（1975）研究了浮力驱动断裂作为岩浆传播的机制.Spence和Sharp（1985）得到了一个无穷大板中流体源作用下流体驱动断裂的相似解，但没有考虑浮力.Spence等（1987）推导出一个弹性固体中半无穷含流体裂纹浮力驱动断裂的稳态解，他们使用了润滑理论描述流体在断裂中的黏性流动，并应用微分方程描述弹性变形.Lister（1990）考虑流体恒常流量下浮力驱动裂纹传播的材料韧性效应，得到了一个裂纹表面位移解，并表明岩脉为泡状头和等宽度长尾巴的形状.Lister和Kerr（1991）描述了岩脉传播中流体和固体耦合效应，指出岩脉中岩浆传输是受浮力和黏性压力降之间的平衡所控制，除了在岩脉尖端外，弹性力起着次级作用.Rubin（1993）得到了一个在给定压力源下岩脉增长的自相似解.Rubin（1995）综述了对围岩破裂包括观察、物理过程和未来方向等各个方面的内容.Meriaux和Jaupart（1998）考虑有限弹性板中浮力和给定超压力驱动下岩脉传播.Bonafede和Rivalta（1999）使用位错模型研究岩脉周围的应力和位移场.Roper和Lister（2005）进一步研究含流体裂纹在超高压半无穷弹性板的传播过程，得到了一个临界长度，即当岩脉长度小于这个临界长度时，岩浆房超高压起着主要作用，当岩脉长度大于这个临界长度时，浮力起着主要作用.Chen等（2007）利用摄动法计算了地壳密度逐渐变化条件下岩脉稳态传播过程，并分析了裂纹稳态传播的条件.

前人有利用弹性半空间的点源模型来模拟研究火山区岩浆囊膨胀与地表变形之间的关系（Mogi，1958），有利用黏弹性围岩球壳模型研究球状岩浆囊周围应力应变响应（Dragoni and Magnanensi，1989），及利用黏弹性扁椭球模型研究岩浆囊膨胀对断层的影响（Newman etal.，2006）.本文所研究的岩浆沿已有断层运移主要涉及黏性岩浆流体在压应力条件下打开已有断层并在弹性地壳中上升的过程，我们采用三维弹性上地壳模型进行模拟.

已有分析表明，在构造应力和岩石的自重作用下，壳内以分异演化后的中酸性岩浆为主的岩浆与围岩之间存在密度差，岩浆会受到一个向上的浮力.由此可能导致上覆围岩中已有断层受到接触岩浆的挤压应力时，断层张开，岩浆向上运移.在此过程中，静岩压力、构造应力以及密度差造成的浮力会随之减小，有可能造成岩浆运移的终止.如果岩浆受到的超压力足够大，则岩浆能够向上运移直至地表.

如上所述，过去关于岩浆在地壳中运移的研究多采用断裂力学方法，研究岩浆驱动下围岩破裂的岩浆传播过程，且大多是二维的.本文针对岩浆在已有断层中运移过程，采用三维有限元法与有限差分法相结合，计算模拟岩浆在地壳中的传播过程.我们利用自行编制的相关程序，对长白山天池火山进行了参数试验的模拟计算.有限差分方法用来分析断层中岩浆的流动，有限元方法用来计算岩脉传播过程中已有断层的打开导致岩体变形过程.

2 物理模型与研究方法

岩浆沿着已知断层面向上运移便形成了岩脉.本文用两个相邻的正方体模型来模拟岩体，相邻的接触面作为断层面，模拟岩浆沿已有断层运移的过程.在水平方向上岩体受到周围构造应力的作用；岩浆在断层中受到构造应力与弹性形变产生的压应力，以及静岩压力与自身重力.在三维模型下，岩浆在断层中运移的宽度（沿断层面）是有限且不变的.

考虑从岩浆房中垂直向上进入围岩的岩浆，假设岩浆为不可压缩的黏性流体，且在运移过程中不会发生固化，利用流体力学的润滑理论描述侵入断层中的岩浆流体（Rice，1968），在流体充满的空间中，根据质量守恒，有

这里w是断层面张开位移，q是流体流量.上式写成差分格式为

其中wi-w0i是第i个断层单元当前张开位移与Δt时间前张开位移之差，Δqi是第i个单元的流量差分.

根据流体流量与孔隙压力的Poiseuille′s定律，有

其中，μ是岩浆的动力黏度，q是岩脉中岩浆流量，w是岩浆侵入使断层张开的位移，p是岩脉中岩浆受到的压应力，ρm是岩浆密度，z表示沿岩脉向上运移方向.驱动岩脉传播压力，即为使岩石弹性形变产生的压应力，包括岩浆流体的压应力、岩石围压和构造应力，因此式（3）中

其中pe＝pe（z）为岩浆压力，使围岩产生弹性形变，ρs是围岩密度，Δσ为底部垂直于岩脉平面的构造应力（Rubin，1995；Roper and Lister，2005）.（4）式可写为

其中

（3）式成为

其中

由（7）式，有

故

同样，

其中pe的导数写为差分格式为

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

结合（2）式，得

其中

对于（14）式，考虑弹性压力pe（z）与张开位移w（z）之间的弹性关系，当给定一组张开位移值时，可以确定一组弹性压力.其中的w初始给定为0.005m，pe两端的条件为

其中S0为侵入岩浆底部构造应力与静岩压力之和，ΔP为给定的超压力，来自于岩浆流体的内部压力与浮力，本研究采用恒常的超压值.S（m）为随运移过程变化的侵入岩浆尖端构造应力与静岩压之和.

对于非连续变形分析中的接触问题，当pe（z）已知，w（z）可根据有限单元方法（王勖成，2003）求出.假设地壳中存在一断层，当侵入岩浆末端的压力pe（a）大于静岩压与构造应力之和时，已存在的断层会张开，岩浆便会继续侵入并向上运移.具体计算时，先由已有断层初始的微小张开位移w0，利用有限差分格式计算出平衡公式（14）中的应力项，得到岩浆压力pe，从而获得围岩的应力状态.得到围岩-岩浆接触面上的边界应力后，用有限单元方法得到此应力边界条件下围岩中各单元的应变量，进而获得接触面上的位移值，即已有断层两侧的张开位移量wi.如果每次获得的接触边界上最上端第m个单元端部节点的张开位移值wm+j大于零，则可以判断断层继续向上张开，否则，岩浆运移停止.如果断层进一步张开，岩浆侵入，则以新获得的wi作为下次计算中的w0，重复上述计算，直到断层面停止张开为止.采用有限单元法计算接触面的位移值，相当于计算岩脉前端岩石所受应力.当其受拉应力作用时断层张开，岩脉向上运移.上述过程将包含已有断层的弹性问题转化为了一个接触问题.

3 计算模型

考虑一个含有已知断层的区域岩体，构建一个三维长方体区域，底面为共用一条边的两个正方形，正方形边长为10km.长方体高h（h＝5000、8000、10000m），构成相互接触的两个正方体区域（图1a）.如图划分网格单元，x方向（断层面方向）中心适当加密，y方向均匀分割，z方向固定为49等份（50层节点），共计3528个单元和5000个单元节点.两正方体的接触面模拟为断层面，在区域下方有一岩浆房，设置初始进入断层微小张开的岩浆位于长方体的底面中心单元并向上延伸5层节点单元.三维模型下岩浆充满断层一个单元的宽度（图1中沿x轴方向）.

图1 （a）计算模型网格，底部红色代表侵入断层的岩浆；（b）侵入岩浆三维示意图，w为岩浆侵入使断层两侧张开的位移；（c）模型受构造应力（Δσ）示意图，p0为底部节点的岩浆压力；（d）侵入岩浆二维示意图，w为断层两侧张开位移，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

模型的底面节点垂直方向固定，四周各面加载10MPa左右的构造应力，且x方向上略大于y方向.构造应力与静岩压力之合力沿垂直方向线性变化.断层从下部开始张开5层节点的初始微小位移w0＝0.005m，最下层处节点的岩浆压力为p0，见式（16）.超压力ΔP为固定大小，当模型深度h一定时，p0恒定不变.

本文所采用的围岩弹性模量E＝50GPa，泊松比ν＝0.25，重力加速度g＝9.8m·s-2.我们对不同组参数（包括岩浆房超压力，岩浆与围岩的密度差及构造应力）组合进行了数值计算.

4 计算结果与讨论

在原有三维模型的基础上，改变设置得到二维模型并进行相关计算，用以对比三维的情形.二维模型采用类似平面应变的模型，仅仅将三维模型中一个单元内的岩浆宽度扩展为沿x方向（断层面走向）上充满断层面，其他条件相似.

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

3-D模型中岩脉的张开位移为3～5m，向上运移逐渐呈现火把形截面；2-D模型张开位移为10～120m，截面形状为盾形或梯形（图2）.两个模型中，岩脉在接近地表时上升速度都会加快，这是因为浅层静岩压力和构造应力均减小之故.2-D模型计算得到的张开位移量与根据断裂力学二维计算的结果一致（Chen etal.，2007）.

表面（地表）垂直位移的形态上，3-D模型先由断层两侧缓慢隆起，进而中心渐渐凹陷，两侧继续突出，使断层与断层两侧突出之间高差增大，最终断层中心较快速突出，形成山峰状；2-D模型则先由断层两侧缓慢抬升，进而中间凹陷成低地，两侧继续抬升，使凹陷与两侧抬升之间高差增大，最终在断层上方较快速地突出，形成山脊状.三维模型中地表最大垂直位移值在0.3m以内；二维模型中地表最大垂直位移值达到12m（脊）和20m（槽）（图3）.二维模型中的地表位移总体上明显大于三维模型中的，且在前期位移量均较小（3-D中为mm～cm级，2-D中为dm～m级），在岩脉接近地表时垂直位移量增大.二维模型中的垂直位移比三维模型的垂直位移大约两个数量级.

两种模型的pe值曲线基本相同（图4），均是先稳定增大，在接近岩脉顶端附近迅速减小至S（m），即构造应力与静岩压力之和（见式（16））.这是由于在每一步岩脉运移之后给定的应力条件限制的.

三维模型中如果限定侵入断层的岩浆宽度适当小时，可以模拟断裂带的火山岩浆运移过程.二维模型可以近似模拟洋中脊岩浆上涌的过程.

我们进行了不同参数下的三维模拟计算，研究了岩浆黏度、密度差、模型深度对最小超压力（岩浆运移到地表所需的最小岩浆超压力）的影响.结果表明，对于深度为10km的模型，相应底部边界的围压（构造应力与静岩压力）为σx＝100MPa，σy＝91.5MPa，岩浆的动力黏度取为2×10-7～2×10-3MPa·s.当超压力达到17～20MPa时，可以驱动岩脉运移到地表.同时，岩浆动力黏度越大，使岩脉运移到地表需要的最小超压力就越大（见表1）.随着模型深度的增加，黏度增大对最小超压力的影响减弱，即增大黏度后，所需最小超压力随模型深度增加缓慢，且在一定黏度范围内，最小超压力变化不大.密度差（围岩有效密度与岩浆密度之差，见式（8））取值为300～700kg·m-3，影响趋势为密度差越小，所需最小超压力越大，且影响幅度比较小——小于1MPa／（0.1g·cm-3）.

对于深度小于10km的模型进行的计算表明，当岩浆囊埋藏较浅时，使岩脉发生运移并运移到地表的最小超压力要比上述10km情形下有明显减小.也许是因为较浅部构造应力和静岩压力的影响减小，使得开始发生运移及维持运移都较为容易.

5 长白山天池火山的应用

长白山天池火山是长白山火山区最具代表性的三座巨型碱性火山之一，是该区发育最为完全的火山之一，是我国境内保存最好的新生代多成因中央式火山（魏海泉，2010）.天池火山从中晚更新世至今有过多次喷发，最近一次大喷发在1215年前后，其规模和程度被认为是近2000年来地球上最大的爆炸式火山喷发之一.自1985年开展天池火山的地震监测以来，2002年7月以前长白山天池火山地区火山地震活动性一直比较低，从2002年7月开始到2004年9月，该地区地震活动性明显增强，地形变也出现了较为明显的变化，这些现象都表明长白山天池火山地区地下岩浆或热液活动正在发生变化（明跃红等，2006），相应于天池火山地下浅层壳内岩浆囊的压力与体积变化.

表1 不同参数下使岩浆运移到地表的最小超压力Table 1 The least overpressure to drive magma up to surface under different parameters

前人对长白山火山区的地质结构进行了深入的研究，对于该地区地壳内岩浆囊的赋存有了进一步认识（刘若新等，1998a；樊祺诚等，2005；张成科等，2002；汤吉等，2001；杨卓欣等，2005；赵金仁等，2003）.近年来的一些研究显示出该地区岩浆囊的新特点.张先康等（2002）利用人工地震测深得到的分层岩浆系统最浅位于8～9km；吴建平等（2007）对天池火山区震群活动的研究指出，5km深度附近可能存在岩浆热液活动或与岩浆增压有关；朱桂芝等（2008）利用Mogi模型联合GPS观测数据并结合遗传算法反演得到长白山火山区岩浆囊位于深度9.2 km处；陈国浒等（2008）基于InSAR、GPS形变场的Mogi模型反演得到的结果推测出天池火山下方7.9km存在一个岩浆压力源；吴建平等（2009）利用面波层析成像得到天池火山附近S波低速层顶层埋深在8～9km；魏海泉（2010）参考有花岗质岩石提出的模型，指出残余的粗面质熔体在5km或更深处的地壳内聚集.这些研究显示出长白山火山区地壳浅部的岩浆囊或许存在于5～10km处.故此我们分别以h＝10km和h＝5km的模型来模拟天池火山的岩浆运移并进行相关讨论.

天池火山岩浆是经过部分矿物结晶分异演化的岩浆，因此其幔源的性质减弱，而偏向于中、酸性的粗面岩质及碱流岩质类型（刘若新等，1998b）.从火山活动的喷出物来看，由早期（上新世到早更新世）的造盾，到后期（早更新世到晚更新世）的造锥，再到全新世-近代以来的浮岩喷发，逐渐由偏基性的玄武质岩石为主过渡到偏中、酸性（粗面岩、碱流岩）的粗面-碱流质岩石为主（刘若新等，1998a；刘若新等，1998b；李春锋等，2006）.从岩石化学组成上看，早期造盾时期以高镁铁性质的拉斑系列岩石为主；后期造锥时期，以碱性、钙碱性的粗面岩、流纹岩为主.浮岩的喷发也以粗面质和碱流质居多（刘若新等，1998a；樊祺诚等，2005）.岩浆的分异演化伴随着各个时期的火山活动.天池火山岩浆也存在着混合作用.来自年代学的研究显示，粗面玄武岩的喷发从早期造盾阶段结束后，一直持续到近代的喷发.对来自该区域的岩石及深浅条带混合等现象的分析表明（樊祺诚等，2005），天池火山造锥喷发的不同时期，存在着粗面质岩浆混合现象.有人认为千年大喷发的触发机制可能是深部的地幔玄武岩浆注入地壳岩浆房的结果（刘若新等，1998b）.水含量的多少对岩浆黏度有显著影响.刘若新等研究发现千年大喷发中碱流质岩浆中的水含量已达到饱和甚至过饱和（刘若新等，1998b），这些因素都会使浅部岩浆的黏度变得更小（Freundt and Schmincke，1992；Hui，2008）.考虑粗面质岩浆的混合以及可能的不同水含量影响，岩浆黏度值可达到104～109Pa·s.

图2 岩浆运移上升图左侧：3-D模型结果，右侧：2-D模型结果；由下至上：岩浆运移随时间变化截图，横轴为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两种情形计算了岩浆发生运移达到地表时所需要的最小超压力，同时模拟了岩浆运移过程中地表随时间的变化特征，并根据岩浆房压力可能由于岩浆流出而导致减小的理论关系，估算了相应情形下使岩浆运移至地表时地下岩浆囊所需的体积大小.

图3 地表垂直位移变化图左侧：3-D模型结果，右侧：2-D模型结果；由下至上：地表垂直位移随时间变化截图，色标为位移值.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），可以看出较浅模型需要的最小超压力要小.前人得到的天池火山千年大喷发时岩浆房超压为6.25MPa和15MPa（魏海泉，2010），前者对应于5km深的碱流质层状岩浆房，后者对应粗面质层状岩浆房（应位于更深处）.较浅部岩浆囊受到的构造应力较小，达到运移条件所需的最小超压力也较小.

地表位移变化与上述模型中给出的结果相似（参考图3-左图）.由地表位移变化图可以看出，地表的隆起随着岩浆向上运移而发生变化，仍然表现出断层面对应岩脉两侧先缓慢隆起，而后向岩脉中心聚拢.地表垂直位移没有呈现中心放射状，原因首先在于考虑到模型的简化以及计算的条件，本文模型设置单一已有断层，而长白山天池火山处于该地区两大主要断裂——北东向的六道沟—天池—甑峰山断裂和北西向的白山镇—天池—金策断裂——的交叉点，这两条走向近乎垂直的断裂控制着天池火山的产出（刘明军等，2004）；其次模型上部没有预设破火山口状的隆起，而对于天池火山已有中心隆起的地形，这类效应或许不会很大.从岩浆运移到地表后的最终形态来看，发生较大垂直位移的范围大致呈现椭圆的形态，也与上述原因有关.

图4 岩浆压力pe变化图左侧：3-D模型结果，右侧：2-D模型结果；由下至上：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.

在本文模型的条件下，排除大规模构造运动使得构造应力发生剧烈变化的情况，若给定适当的超压力，则可以粗略估算这种运移发生时天池火山下方浅部岩浆囊的大小.可以利用下列公式（Maccaferri etal.，2010）估算模型下方岩浆囊的容量变化与超压力变化的关系，

其中V为岩浆囊中的岩浆体积V0与侵入岩脉中的岩浆体积Vout之和，Kp为岩浆的体积模量，ΔP为超压力.上式可改写为

根据上式将模型中的超压力设置为随着岩浆向上运移而发生变化.对于不同深度的模型取最小超压力作为初始ΔP0，即5km时ΔP0＝10MPa；10km时ΔP0＝21MPa.岩浆体积模量Kp＝10000MPa.

由每一步岩脉向上运移的高度可算出该次运移由岩浆房向已有断层中侵入的岩浆体积，即为Vout.为保证岩浆可以在给定初始超压力下能够运移至地表，计算出应给定的在运移发生前岩浆囊的体积.对于5km的模型，如果岩浆囊超压力为10MPa，下方岩浆囊体积大于6×109m3时可以保证岩浆运移到地表；对于10km的模型，如果岩浆囊超压力为21MPa，下方岩浆囊容量为2×1010m3时可以保证岩浆运移到地表.如果考虑层状岩浆囊，假设半径为2～6km，则此时5km模型岩浆囊的厚度约为50～500m，10km模型岩浆囊厚度约为200～1600m.

上述对岩浆囊大小的估计，是在构造应力不发生剧烈变化且给定超压力情况下所作的计算，如果岩浆能够运移到地表，那么天池火山下部对应深度可能存在的岩浆囊原有的大小.

6 结论

对于10km深度的岩浆囊，考虑岩浆浮力和构造应力情况下，岩浆黏度在2×10-3到2×10-7MPa·s范围时（相当于玄武岩浆），当超压力达到17～20MPa时，岩浆可以沿着已有断层向上运移至地表，从而形成岩脉.在总体上看，岩浆黏度越大，岩脉发生运移且至地表的最小超压力就越大；密度差（围岩有效密度和岩浆密度之差）越小，所需最小超压力就越大.

岩浆囊埋藏越浅（10～8km），使岩脉发生运移且至地表的最小超压力就越小.

三维模型中岩脉的张开位移为3～5m，二维模型中则为10～120m.在接近地表时，两种模型中岩脉的运移速度都会加快.

三维模型得到的地表垂直位移，由中心凹陷、两侧突出逐渐形成中心突出的峰状；二维模型中则由中间凹陷为槽状逐渐形成中间凸起的脊状.三维模型可以近似模拟断裂带控制的火山岩浆活动，二维模型可以近似模拟洋中脊的情况.

用本文模型模拟长白山天池火山岩浆运移时，地表垂直位移没有呈现中心放射状，原因在于设置的单一已有断层，而天池火山处于该地区两大断裂的交叉点，这两条走向近乎垂直的断裂控制着天池火山的产出；其次模型上部没有预设破火山口状的隆起，而对于天池火山已有中心隆起的地形，这类效应或许不会很大.

若假设长白山天池火山下5km或者10km深度存在岩浆囊，应用本文的模型，如果岩浆要运移至地表，超压力至少应为10MPa和21MPa.当构造应力不发生剧烈变化且给定超压力分别为10MPa和21MPa时，若岩浆运移到地表，岩浆囊应有的最小体积约为6×109m3（5km深）和2×1010m3（10km深）.相对应于半径为2～6km的层状岩浆囊时，厚度约为50～500m（5km深）和200～1600m（10km深）.

本文对于长白山天池火山地区岩浆活动机制和特征及喷发条件的研究结果为未来该地区火山喷发预测和监控提供了参考.

致谢 感谢中国科学院地球深部研究重点实验室提供计算模拟平台.感谢中国科学院超级计算中心提供科学计算网格服务.

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，Fan 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，Fang 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形变场的长白山地区火山岩浆囊参数模拟研究.地球物理学报，51（4）：1085-1092，doi：10.3321／j.issn：0001-5733.2008.04.017.

樊祺诚，隋建立，孙谦等.2005.天池火山千年大喷发的岩浆混合作用与喷发机制初步探讨.岩石学报，21（6）：1703-1708，doi：10.3321／j.issn：1000-0569.2005.06.017.

李春锋，张兴科，张旸等.2006.长白山天池火山的地质构造背景.地震地磁观测与研究，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.长白山天池火山近代喷发.北京：科学出版社.

明跃红，苏伟，房立华.2006.长白山天池火山地震类型及火山活动性的初步研究.中国地震，22（1）：56～63，doi：10.3969／j.issn.1001-4683.2006.01.006.

汤吉，邓前辉，赵国泽等.2001.长白山天池火山区电性结构和岩浆系统.地震地质，23（2）：191-200，doi：10.3969／j.issn.0253-4967.2001.02.008.

王勖成.2003.有限单元法.北京：清华大学出版社.

魏海泉.2010.长白山火山岩浆柱岩浆上升作用过程.地学前缘，17（1）：11-23.

吴建平，明跃红，张恒荣等.2007.长白山天池火山区的震群活动研究.地球物理学报，50（4）：1089-1096，doi：10.3321／j.issn：0001-5733.2007.04.016.

吴建平，明跃红，苏伟等.2009.长白山火山区壳幔S波速度结构研究.地震地质，31（4）：584-597，doi：10.3969／j.issn.0253-4967.2009.04.002.

杨卓欣，张先康，赵金仁等.2005.长白山天池火山区三维地壳结构层析成像.地球物理学报，48（1）：107-115，doi：10.3321／j.issn：0001-5733.2005.01.016.

张成科，张先康，赵金仁等.2002.长白山天池火山区及邻近地区壳幔结构探测研究.地球物理学报，45（6）：812-820，doi：10.3321／j.issn：0001-5733.2002.06.008.

张先康，张成科，赵金仁等.2002.长白山天池火山区岩浆系统深部结构的深地震测深研究.地震学报，15（2）：135-143，doi：10.3321／j.issn：0253-3782.2002.02.003.

赵金仁.2003.长白山天池火山区上地壳三维速度层析成像.地球物理学报，46（6）：796-802，doi：10.3321／j.issn：0001-5733.2003.06.011.

朱桂芝，王庆良，石耀霖等.2008.各向同性膨胀点源模拟长白山火山区岩浆囊压力变形源.地球物理学报，51（1）：108-115，doi：10.3321／j.issn：0001-5733.2008.01.014.