火山活动是地球深部岩浆上涌至地表的一种地质现象，是地球内部热能爆发性释放的一种形式。一座火山从最初岩浆房的形成到后来岩浆的分异、迁移、喷发、固结等各个环节无不与温度密切相关。高温岩浆的产生、赋存和运移一方面受控于特定的岩石圈热构造背景，另一方面必然对岩石圈温度分布产生深刻的影响。因此火山地区，特别是新生代火山活动地区岩石圈温度的时间和空间变化蕴含着丰富的地球科学信息，是了解岩石圈的动力学过程、岩石圈演化、火山活动、地球内部物质和能量迁移、多种矿产资源的形成以及地热资源潜力的重要窗口。

火山区岩石圈温度场有如下特点：(1) 无高温异常体的岩石圈温度场主要取决于稳态地温场，稳态地温场主要是由地球内热、岩石放射性热源，地层导热性质及地表温度决定的，等温面大致呈水平分布，温度随深度变化率大致为3℃/100m左右，如果没有构造变动，可在数百万年甚至更长时间大体不变(汪集暘, 2015)。而在火山活动构造区，由于高温岩浆侵入体的存在，导致其温度场除了受稳态地温场控制外，还受到来自高温岩浆所携带的热量的影响(张旗等, 2014)；(2) 有别于稳态地温场，岩浆侵入体的影响范围和持续时间非常有限，唐晓音等(2013)对琼东南盆地长昌凹陷火成岩侵入体的研究表明，一个面积为300km²，厚约10km的侵入体，影响范围约在3km左右，10Ma后侵入体影响消失，温度场恢复至未侵入前的状态。Daniels et al. (2014)对岩墙侵入体的研究表明，平均面积为16km²，厚约10km的侵入体，对围岩的影响时间约在1Ma，影响范围也约在3km左右。

在诸多的岩石圈温度场分布研究方法中，传统勘探方法成本较高，且受制于很多因素。而随着计算机技术应用而快速发展起来的数值模拟方法，成本相对较低，更为重要的是在结合了地质，地球物理，地球化学，火山学等学科的成果，采用合理的模型和计算方法，可以解释火山区岩石圈的地温场的演化，从而弥补实际观测和实验数据的不足。国外已经有很多学者(Carslaw and Jaeger, 1959; Crank and Nicolson, 1996; Douglas and Rachford, 1956; England and Thompson, 1984; Gartling and Thomas, 1984; Finsterle et al., 2014; Furlong and Chapman, 1987; Mundry, 1970; Rutqvist et al., 2002; Wang and Hyndman, 1993)在地热数值模拟基础理论和案例实际模拟计算分析方面做了大量的工作。国内学者(熊亮萍和张菊明, 1988; 唐晓音等, 2013; 王民等, 2010; He, 2015; Sun et al., 2013)也在岩石圈温度场分布数值模拟方面取得了一些初步成果，但对于火山区岩石圈热结构的数值模拟研究案例还比较少。

考虑到火山地区地热地球物理建模、计算流程、参数设定的复杂性，本文利用在工程热物理领域广泛应用的有限元软件ANSYS WORKBENCH进行建模、计算和分析。相较于其他计算软件，ANSYS WORKBENCH在地球科学模拟方面具有如下优势：1) 采用有限元方法，适用于地质模型中多不规则形状的网格化处理；2) 丰富的建模接口(Hatch, 2000)，可以和比较通用的CAD建模工具如SOLIDWORKS、Pro/e、AutoCAD (Stolarski et al., 2011)等进行连接；3) 建模流程化，对于火山区多期次喷发问题可以加入多个子过程，使模型更加接近实际情况；4) 可以设定多组非线性参数，降低结果误差；5) 集成多个高精度的计算模块，如Multiphysics、flunnt、CFX等，可以求解复杂的物理过程；6) 可加入APDL语言(陆旭等, 2009)来灵活的处理模型，控制计算，并且可对计算结果进行多维度提取。

本文调整了ANSYS WORKBENCH的各项设置，加入多组APDL语言控制计算流程使其适用于火山区岩石圈温度场研究，并建立一系列简单的柱状模型，模拟了不同规模、不同顶板埋深的侵入体热扩散过程，并将上述模型用于长白山天池火山千年喷发热扩散过程的模拟，构建了上中地壳岩浆房和上地壳岩浆房两个不同情景的长白山天池火山岩浆侵入概念模型，模拟了在这两种情景下的岩石圈温度场分布变化。虽然本文在模拟计算中对岩浆体的侵入和冷却过程作了高度简化，但这将为今后把ANSYS WORKBENCH应用于复杂条件的火山区岩石圈温度分布研究奠定基础。

笛卡尔坐标系下三维稳态导热方程为：

(1)

式中λ_i(i=x, y, z)为导热系数，W为源汇项，T为温度，在本计算中W为岩石放射性生热率，单位为W/m³。假设岩石导热性能为各向同性的，则上式简化为：

(2)

相应的非稳态导热方程简化为下式：

(3)

式中t为时间，ρ为密度，c为比热容。

ANSYS采用有限元方法对稳态传热方程进行离散，可以得到n个联立的线性代数方程组，经过整理，可以得到下式：

(4)

式中[K]为n维热传导矩阵，包含导热系数、对流系数、辐射率和形状系数；[T]和[Q]分别为n个节点的温度和热流向量。

瞬态传热过程是指一个系统的加热或冷却过程。在这个过程中系统的温度、热流率、热边界条件及系统内能随时间都有明显变化。根据能量守恒原理，瞬态热平衡可以表达为：

(5)

式中[C(T)]为随温度变化的比热矩阵，[T′]为温度对时间的导数，[K(T)]为随时间变化的传导矩阵，[T]为温度列阵，[Q(T)]为节点随温度变化的热流率，包含热生成。

本文计算流程大体分为两个阶段：第一阶段模拟在无高温岩浆体侵入情况下的稳态地温场分布，第二阶段以稳态地温场分布为初始条件，计算分析高温岩浆侵入后岩石圈温度场的时空分布。

本文柱状侵入体的地质模型为直径40km、深度40km的圆柱体，假设侵入体位于模型中心(图 1中绿色部分)，正上方另设一底部半径1.2km、顶部半径1km、高1km的圆台以模拟火山锥。坐标参考平面为地面。所建模型依据侵入体半径0.5km及1.0km分为两组，每组分别又有火山颈(岩浆喷出火山锥)及顶板埋深分别为0、2、5(岩浆未喷出)4个子系列，总共8个子模型。由于其对称性，本文只取模型的1/8作为研究对象(图 1)。

图 1

Fig. 1

图 1 模型的几何形态及边界条件设定 Fig. 1 The model geometry and the boundary condition

边界条件设定如下：模型下边界为恒热流边界，取地幔热流值为40mW/m²。上边界地面为恒定温度边界，基准面地表温度为10℃，火山锥侧面边界为恒定温度边界，温度随高程的变化率为5℃/km。火山锥顶岩浆侵入体与空气接触面为对流边界，岩浆和空气的换热量Φ可由下式得到，

(6)

式中A为接触面积，t₁为岩浆侵入体上表面温度，t₀为空气参考温度，设为5℃，h_m为不同岩浆温度下的平均表面传热系数，可由下式得到(Patankar, 1980)

(7)

式中u_∞为空气速度，取为10m/s，v为空气的运动粘度，Pr为空气的普朗特数，λ为空气热导率，v，Pr，λ随岩浆温度的变化值可查找干空气性质表得到(Bejan, 2013)，图 2a是模拟计算中采用的火山锥顶面平均表面传热系数。

图 2

Fig. 2

图 2 火山锥顶面平均表面传热系数随岩浆侵入体上表面温度的变化(a)和岩浆侵入体密度随温度的变化(b) Fig. 2 Upper surface temperature of the magma chamber versus the average surface heat transfer coefficient at the top of the volcanic cone (a) and temperature versus density of the magma column (b)

在计算模型的热物性参数取值方面，综合参考Hasterok and Chapman (2011)所做的工作及Wang and Cheng (2012)计算中所采用的数据，取地壳厚度为40km，分为上、中、下三层，厚度分别为15km、15km和10km，采用的热物性及生热率参数如表 1所示。计算过程中岩浆侵入体的密度随温度是非线性变化的，如图 2b所示。

表 1

Table 1

表 1(Table 1)

表 1 本文采用的岩石热物性参数(据Hasterok and Chapman, 2011; Wang and Cheng, 2012) Table 1 Thermal physical properties of rocks used in this study (after Hasterok and Chapman, 2011; Wang and Cheng, 2012) 地层 密度(kg·m-3) 热导率(W·m-1· K-1) 比热(J·kg-1· K-1) 生热率(μW·m-3) 上地壳(0~15km) 2750 2.6 1050 1.5 中地壳(15~30km) 2900 2.5 970 0.4 下地壳(30~40km) 3200 2.8 950 0.02 岩浆侵入体 变化(图 2) 2.5 1100 0.002

表 1 本文采用的岩石热物性参数(据Hasterok and Chapman, 2011; Wang and Cheng, 2012) Table 1 Thermal physical properties of rocks used in this study (after Hasterok and Chapman, 2011; Wang and Cheng, 2012)

由于岩浆侵入体最终冷却温度会趋同于同深度稳态地温场温度，而稳态地温场的温度是深度的函数，因而在冷却过程中侵入体与围岩的温度趋于平衡的程度随深度而变化。Mundry (1970)在研究岩浆冷却过程提出了利用无量纲的冷却温度来描述侵入体的冷却程度。其中T_稳态为稳态温度，T_侵入为岩浆侵入初始温度，T为冷却过程温度。本文在此基础上引入冷却指数α

(8)

即岩浆侵入体内某一位置冷却过程温度与稳态温度之差与初始侵入温度之差的比值来描述侵入体的冷却程度，若冷却指数大于95%即认为冷却过程基本结束。

由于柱状侵入体的传热方式是由中心向外扩散的，温度随着与对称轴的水平距离(径向距离)的增大而降低，因而在确定侵入体某一深度的冷却时间时，以距离为0处，即中心线处冷却指数大于95%的时间为准。如侵入体给定深度位置的稳态温度T₀为50℃，若侵入体的初始温度T_侵入=1200℃，则当岩浆侵入体的中心位置温度降低到107.5℃时便可认为这一深度的岩浆侵入体已经基本冷却。

围岩受侵入体热扩散扰动以干扰温差(受扰动后的温度与对应处侵入体未侵入时温度之差)ΔT为依据，当ΔT < 5℃时便认为此处温度场基本不受侵入体热量扩散的影响。

不同规模，不同埋深的岩浆侵入体在不同深度处冷却指数随时间的变化如图 3所示。为较直观的获得冷却时间，本文对图 3中冷却指数大于95%的部分做了放大处理。由图 3分析可知，深度增加，冷却时间随之增加；半径由0.5km增加到1km，冷却时间显著增加。表 2是图 3所示8个子模型冷却时间的汇总。结果表明，岩浆的冷却时间与顶板埋深成负相关关系，与侵入体的规模(半径)平方约成正比关系，与测点深度成正相关关系，但到一定深度趋于定值。例如，对于半径为0.5km、顶板埋深为0km的侵入体，深度从0.5km增加到7km，冷却时间从23.0ka增加到50.6ka但从7km以后，冷却时间都为50.6ka；而0.5km半径和1.0km半径系列达到基本冷却的时间分别约为50ka和200ka，冷却时间大致与半径平方成正比关系。

图 3

Fig. 3

图 3 处于不同深度不同规模和埋深的柱状岩浆侵入体冷却指数随时间的变化 Fig. 3 Variation of cooling factor with time at given depths for different scales and intrusion depths of a magma column

表 2

Table 2

表 2(Table 2)

表 2 火山活动侵入岩浆达到冷却指数为95%所需要的时间(单位ka) Table 2 Lapse time (ka) for the cooling factor to be equal or greater than 95% for different scales and intrusion depths of a magma column at given depths 深度(km) 岩浆侵入柱体顶板埋深(km) 侵入体半径0.5km 侵入体半径1.0km 火山颈 0 2 5 火山颈 0 2 5 0.5 28.0 23.0 - - 74.2 58.0 - - 1.0 32.7 28.4 - - 92.5 78.4 - - 1.5 42.7 33.5 - - 116.7 94.2 - - 2.0 43.4 34.0 21.2 - 131.2 107.5 86.2 - 4.0 50.4 49.8 37.2 - 182.5 152.1 135.4 - 5.0 50.8 50.4 42.9 21.9 200.6 185.2 153.5 85.6 7.0 50.8 50.6 50.4 50.4 201.5 198.4 180.2 137.3 10.0 50.8 50.6 50.4 50.4 201.5 198.4 197.9 194.2

表 2 火山活动侵入岩浆达到冷却指数为95%所需要的时间(单位ka) Table 2 Lapse time (ka) for the cooling factor to be equal or greater than 95% for different scales and intrusion depths of a magma column at given depths

侵入体对围岩温度场影响如图 4所示。图中左侧白线所框部分为侵入体所在位置，且本文认为火山锥中包围火山颈部分也为围岩的一部分。高温岩浆侵入以后，侵入体逐渐把能量传递给围岩，温度不断降低，传热温差会不断缩小，最终冷却温度会趋同于同深度地温场温度。根据模拟计算结果，可以看出以下特征：

图 4

Fig. 4

图 4 侵入体侵入后不同模型不同时刻温度场剖面 Fig. 4 Cross sections of temperature distribution for different models at given lapse times after intrusion

(1) 对于火山颈半径为1km侵入体直达地表的模型，从岩浆柱侵入开始(图 4a)到0.5Ma (图 4c)，侵入体温度迅速下降，以侵入体为中心，围岩等温面逐渐向上抬升，侵入体对围岩的横向热扰动范围逐渐从侵入体边界扩展至7km，纵向范围从火山锥顶逐渐减小至2km深度，1.4Ma后(图 4e)，侵入体对围岩的热扰动基本结束，侵入体与围岩温度趋于一致，温度场逐渐恢复至稳态地温场。

(2) 半径1km顶板埋深5km的柱状侵入模型计算结果(图 4f-i)表明，从侵入开始到0.1Ma，侵入体温度迅速下降，侵入体对围岩的横向热扰动范围逐渐增加至5.5km，纵向热扰动范围从顶板5km深度处逐渐扩展至1.7km深度；从0.1Ma到0.5Ma，侵入体温度继续下降，横向热扰动范围从5.5km增至7km，纵向热扰动范围由1.7km深度回落至3.5km深度；1Ma后，侵入体对围岩的热扰动基本结束，侵入体与围岩温度趋于一致，温度场逐渐恢复至稳态地温场。

(3) 对比火山颈半径为0.5km的较小规模岩浆侵入冷却过程模拟模拟结果(图 4j-l)可以看到，从侵入开始到0.1Ma，侵入体温度迅速下降，岩浆侵入体对围岩的横向热扰动范围从侵入体边界逐渐增加至3km，纵向范围从火山锥顶逐渐减小至0.7km深度，0.5Ma后，侵入体对围岩的热扰动基本结束，在同一深度侵入体内外温度趋于一致，地温场逐渐恢复至稳态。

(4) 对比侵入体半径同为1km，顶板出露地表(火山颈1km)和顶板埋深为5km的两个模型计算结果(见图 4a-e与图 4f-i)可以发现，埋深大小的影响主要体现在侵入体对围岩的纵向热扰动范围及持续时间上，而对横向热扰动规律影响很小。若岩浆喷出，随着时间的推移，侵入体对围岩的纵向热扰动范围会从火山锥顶逐渐减小；若岩浆未喷出，侵入体顶板埋深一定时，侵入体对围岩的纵向热扰动范围会在一段时间内从顶板处增加至一定深度，而后再逐渐减小直至消失。

(5) 对比半径1km火山颈与半径0.5km火山颈两组模型计算结果(见图 4a-e和图 4j-l)，侵入体规模决定了其对围岩热扰动的范围及持续时间。比如1km火山颈模型中侵入体对围岩温度场影响可持续至1.4Ma，最大横向热扰动范围可至7km，而在0.5km火山颈模型中，热干扰持续时间只有不到0.5Ma，最大横向热扰动范围仅为3km。侵入后0.1Ma，1km火山颈模型的侵入体对围岩的纵向热扰动范围只是火山锥体基本冷却。而在0.5km火山颈模型中，火山喷发0.1Ma以后，侵入体对深度0.7km以上的浅层地层的热干扰已经基本消失。

(6) 从火山颈半径为1km的模型计算结果(见图 4a-e)可以看到，从地表到某一深度，侵入体对围岩的热扰动随着深度增大影响范围逐渐增大，达到最大值后随着深度增加，影响范围逐渐减小。产生这个现象的原因在于上表面为恒温边界，且地表温度远远低于侵入体及受扰动的高温围岩，因此越靠近上表面，散热越剧烈。这个现象与实际情况是相符的，由于地面是岩石圈和大气圈的热交换边界，浅层地表温度会受到大气与地下温度场交互影响而基本维持在一个稳定范围内(Huang et al., 2000)，且远远小于侵入体及受扰动的高温围岩，因此根据传热理论，在传热能力基本不变的情况下，越接近地表，传热温差越大，侵入体传入围岩的温度就越容易散失到空气中，从而对近地表围岩的影响范围相对较小。

我国东北地区是东亚最重要的现代火山活动区之一，多次发生了以碱性玄武质、粗面质、流纹质为特征的火山喷发，形成大片熔岩台地和一系列新生代火山锥。其中长白山火山乃是我国大陆地区最为活跃、最具潜在喷发危险的火山(刘嘉麒等, 2015)。

火山岩石学证据表明(Liu et al., 2015; Zhang et al., 2015; Zou et al., 2014)，长白山天池火山经历了早更新世造盾、中-晚更新世造锥和全新世喷发三个发展阶段，初始岩浆源自地幔，但喷发前在地壳岩浆房中多有驻留。火山活动具有多旋回、多期次、多阶段喷发、在时间和空间分布上此起彼伏的特征，在全新世有三个喷发期，第一喷发期距今约5000~4000年，第二喷发期为公元946年左右，而第三喷发期距今仅约300年。其中以第二次喷发规模最大，即所谓“千年大喷发”，是全球近2000年来最大的一次喷发事件，喷发的火山灰不仅沉积到我国东北、朝鲜半岛、日本海及日本北部, 甚至在格陵兰冰川中都能检测到(Tian et al., 2011; 刘祥, 2006)，喷发的物质总量达100km³，在主峰山顶形成了直径约5km的盆状破火山口，天池便是火山口积水形成的湖泊。

关于长白山天池火山岩浆房的深度和规模，不同的学者根据不同的资料和研究方法得到的认识有所不同。最近，郭文峰等(2015)采用岩石热力学模拟方法，研究了长白山天池火山岩浆房的成分和埋深，认为长白山火山全新世岩浆来自中地壳上部(上中地壳)的岩浆房，深度大致为15~18km (0.5~0.6GPa)。而利用大地电磁测深数据反演得到的岩石圈电性结构则显示，在天池火山口下方存在明显的直立型岩浆通道，岩浆通道在下方约5~8km位置形成关闭，并与埋深约为7km、电阻率小于10Ω·m的明显低阻异常体对接(仇根根等, 2014)，这一低阻异常体的埋藏深度与根据高精度卫星遥感地形形变数据反演得到的上地壳岩浆房的埋藏深度基本一致(何平等, 2015)。

根据上述研究成果，本文构建了岩浆房分别位于上地壳和上中地壳的两个长白山天池火山岩浆侵入概念模型，岩浆房的形态以一直径为20km、高度为5km的扁球体近似，岩浆通道则以一通过岩浆房中心、贯穿整个模型、直径为5km的垂直柱状体近似。在上地壳岩浆房情景中，假设岩浆房中心距模型基准面深度为10km如图 5a所示；在上中地壳岩浆房情景中，岩浆房中心距模型基准面深度则设定为20km，如图 5b所示。假设岩浆房和岩浆通道的侵入体是瞬时形成的。其它物性参数及边界条件与前面的简单模拟计算相同。

图 5

Fig. 5

图 5 岩浆房分别位于上地壳(a)和上中地壳(b)的长白山天池火山岩浆侵入概念模型 Fig. 5 Geometry of the two magma intrusion models with the magma chamber located in the upper crust (a) and upper-mid crust (b), respectively, for the Changbaishan Tianchi Volcano

图 6是两个模型情景下一千年前长白山天池火山即将喷发的初始(0ka)、现今(1ka)、和11ka及101ka后的地温分布。

图 6

Fig. 6

图 6 长白山上地壳岩浆房及上中地壳岩浆房模型侵入体侵入后不同时刻温度场剖面 Fig. 6 Snapshots of modeled temperature cross sections after magma intrusion for the two conceptual models for the Changbaishan Tianchi Volcano (a) and (e): 0ka, or initial state; (b) and (f): 1ka after eruption, or current; (c) and (g): 11ka after eruption, or 10ka from now; (d) and (h): 101ka after eruption, or 100ka from now; (i), (j), (k), and (l) are zoom-in views of (b), (c), (f), and (g) respectively

从千年大喷发至今，岩浆侵入体温度并无明显下降，只是在距侵入体与围岩边界约250m范围内的侵入体温度降低，最高降幅超过500℃，而侵入体对围岩温度场的扰动范围也仅仅局限在距边界约500m的区域内，紧靠近侵入体的围岩温度增幅最大超过500℃，围岩及侵入体等温线以侵入体与围岩边界为约束发生了变形，且在此范围内形成一个较大的横向地温梯度带。对比两个组模拟计算结果，上中地壳岩浆房模型较上地壳岩浆房模型，侵入体温度冷却程度更慢，其温度降幅最高不超过400℃；侵入体对围岩温度场的扰动也更小，围岩温度最高增幅不超过400℃。

模拟结果显示，距今10ka以后，侵入体冷却范围将扩大到1.1km，上地壳和上中地壳岩浆房两个模型侵入体最高降温幅分别约为600℃和500℃，侵入体对围岩的扰动范围均增至距边界1.9km，围岩温度最大增幅分别为500℃和400℃左右。

距今100ka以后，不论是上地壳岩浆房模型或者是上中地壳岩浆房模型，温度场的变化较10ka时的情景将会有很大的变化，围岩及侵入体等温线的变化将进一步脱离侵入体与围岩边界的约束，岩浆房及邻近区域的等温线将从椭圆向圆形发展。在上地壳岩浆房模型中，如图 6d所示，岩浆房及柱状侵入体内高温800℃以上部分(图中黄色及金黄色所标记区域)明显小于上中地壳岩浆房模型，如图 6h所示，表明上地壳岩浆房模型中岩浆房温度下降速度较上中地壳模型快。侵入体对围岩温度场的扰动范围均为距边界5km区域，且岩浆房附近距岩浆房距离相同处，上地壳岩浆房模型围岩温度增幅均大于上中地壳岩浆房模型，表明虽然上地壳岩浆房模型侵入体温度下降速度较上中地壳岩浆房模型侵入体快，但这并未使其对围岩温度场的扰动范围增长速度发生改变，而是使得上地壳岩浆房模型岩浆房附近围岩温度增幅大于上中地壳岩浆房模型。

本文首先计算了一系列简单的柱状侵入体的热扩散过程作为利用ANSYS模拟火山区岩石圈温度场分布的基础，结果表明，半径1km的柱状侵入体冷却时间约为200ka，半径500m为50ka，冷却时间与半径平方成正比，这与采用数学解析方法所得结果相吻合(Mundry, 1970)。而在侵入体对围岩影响方面，Daniels et al. (2014)计算表明，平均面积为16km²，高为10km的侵入体对围岩的热干扰不超过1Ma，其所研究的侵入体体积为160km³，与本文的1km半径柱状侵入体的105~120km³基本在一个数量级上，所得侵入体对围岩的热干扰不超过1Ma结果与本文的1.4Ma也是吻合的。唐晓音等(2013)计算结果表明，面积300m²，高约10km的侵入体对围岩的热干扰持续时间为10Ma，将本文结果与其对照，本文按其标准，1km半径柱状侵入体冷却时间应约在0.5Ma左右。其研究的侵入体体积为3000km³，约是本文1km柱状侵入体的105~120km³的25倍，若简单的按侵入体对围岩的热干扰时间与体积成正比来计算，本文1km半径柱状侵入体冷却时间应约在0.5Ma左右，因而唐晓音等人的计算结果与本文也是相符的。

上述对比表明，虽然简单柱状模型将火山侵入体模型大大简化，但是仍然能够较好的描述侵入体热扩散过程，因而将其作为利用ANSYS模拟火山区岩石圈温度场分布的基础是可行的，将柱状模型的结构依据具体地区实际情况进行重建，且对计算参数，计算流程，针对具体问题进行优化，并加入新的模拟因素，可为人们了解岩石圈的动力学过程、岩石圈演化、火山活动、地球内部物质和能量迁移、多种矿产资源的形成以及地热资源潜力提供一个新的数值方法。

本文在柱状侵入体模型的基础上进行了模型重建，综合长白山最新研究成果，构建了上中地壳岩浆房和上地壳岩浆房两个不同情景的长白山天池火山岩浆侵入概念模型，以模拟长白山大喷发对于地温场的影响。

两个长白山天池火山模型的模拟结果表明，不论千年大喷发的岩浆房是位于上地壳还是上中地壳，侵入体对现今地温场的影响主要局限在距侵入体和围岩边界至多500m的范围内，在此区域产生了短距离但梯度较大的横向地温梯度，随着时间的推移，受到岩浆侵入体热扩散影响的区域将会逐渐增大，横向地温梯度也会发生相应变化。图 7a和图 7b分别为上地壳岩浆房和上中地壳岩浆房两个模型中围岩5km深度处温度随距侵入体与围岩边界距离变化。如果以围岩受干扰温差ΔT>100℃的区域为高温地热资源富集区，现今高温地热资源主要集中在侵入体及其周边0.3km的范围内，1ka后将增加至0.45km的范围，10ka后至0.85km，100ka后在上地壳岩浆房情景中增至3.3km，在上中地壳岩浆房情景中增至2.1km。

图 7

Fig. 7

图 7 两个模型中围岩5km深度处温度随距侵入体与围岩边界水平距离变化图 (a) 上地壳岩浆房模型；(b) 上中地壳岩浆房模型 Fig. 7 Temperature of the country rocks versus distance from the edge of the intrusion at the depth of 5km (a) the upper crust magma camber scenario; (b) upper-middle crust magma chamber scenario

目前，国内对于火山区岩石圈温度场的研究还比较欠缺，本文利用在工程热物理领域广泛应用的ANSYS WORKBENCH软件，构建了一系列柱状模型来模拟火山岩浆侵入体的热扩散过程，并综合长白山最新研究成果，构建了上中地壳岩浆房和上地壳岩浆房两个不同情景的长白山天池火山岩浆侵入概念模型，模拟了长白山大喷发对于温度场的影响，得出以下几点认识：

(1) 半径1km的柱状侵入体冷却时间约为200ka，500m为50ka，冷却时间大致与半径平方成正比；

(2) 侵入体对围岩的热扰动的范围和持续时间很大程度上取决于侵入体的规模，比如半径为1km的火山颈对围岩温度场的影响可持续1.4Ma，最大横向热扰动范围可至7km，而半径0.5km的火山颈热扰动持续时间只有不到0.5Ma，最大横向热扰动范围仅为3km；

(3) 对长白山天池火山千年大喷发的模拟结果表明，天池火山颈中心的温度目前基本上仍然保持原有的高温状态，温度仅在厚度约250m的侵入体外圈有明显降低，10ka以后冷却范围可能扩展到1.1km，100ka明显的冷却范围将扩展至侵入体最内部。

(4) 不论其岩浆房的位置是在上地壳或是中地壳上部，长白山天池火山千年大喷发对于现今温度场的扰动目前尚局限在距侵入体边界500m的范围内，区内地温变化范围基本上没有脱离侵入体的轮廓线的约束，这可能是火山区内地热显示不强烈的一个重要原因；

(5) 随着时间的推移，与长白山天池火山千年大喷发相关联的岩浆侵入体对区内地温场的影响将逐步得到进一步显现，10ka后侵入体对围岩的扰动范围将增至距边界1.9km，且将在影响区域内形成了一个较大的横向地温梯度带；

(6) 在10ka时间内，上地壳岩浆房模型和上中地壳岩浆房模型两个模型中侵入体对于浅部地温场对影响没有明显差异，在距今10ka以后到100ka这段时间内，不论是上地壳岩浆房模型或者是上中地壳岩浆房模型，区域地温场将发生很大的变化，但在此过程中，上地壳岩浆房模型中侵入体的温度下降速度明显快于上中地壳岩浆房模型；

(7) 火山活动是重要的岩石区构造热事件，火山地区又往往是最具地热开发前景的地区(黄少鹏, 2014)，深化对于火山岩浆活动与地温场的影响的认识，不仅具有重要的理论研究意义，而且对于火山地区的地热资源勘探、开发具有重要的指导意义。