-
花岗质岩浆侵入不仅引起围岩的接触变质作用和变形（Brown and Solar，1999; Hobbs and Ord，2010; Brown，2013），还提供了成矿流体（Wang Guoguang et al.，2020）。岩浆侵入过程中形成的断裂和微裂隙是熔/流体活动和元素迁移的重要通道，也是成矿元素富集沉淀的主要场所。花岗岩与许多金属矿床的形成密切相关，前人围绕花岗岩的构造环境、成因机制和成矿作用开展了大量研究，例如:安第斯造山带和特提斯造山带的斑岩型铜-金矿床（Camus and Dilles，2001; Hou Zengqian et al.，2003; Richards，2015）、冈底斯的斑岩型钼矿床（Zhao Junxing et al.，2012）、秦岭的斑岩型和矽卡岩型钼矿（Chen Yanjing，2010; Yang Zongfeng et al.，2011）、华南的钨锡矿和稀有金属矿（Hu Ruizhong and Zhou Meifu，2012; Mao Jingwen et al.，2012; Chen Jun et al.，2013）、松潘甘孜造山带的锂矿（Fu Xiaofang et al.，2014; Wang He et al.，2017; Xu Zhiqin et al.，2018）、喜马拉雅造山带的稀有金属矿（Wang Rucheng et al.，2017; Wu Fuyuan et al.，2021; Qin Kezhang et al.，2021a）等。
-
此外，岩浆侵位时的热烘烤效应可促进富有机质围岩的生烃作用（Galushkin，1997; Chen Zhenyan et al.，1999; Othman et al.，2001），而接触变质带附近的高温和脆性变形又可能破坏油气资源的保存（Svensen et al.，2007; Agirrezabala et al.，2014; Li Xia et al.，2016; Zhang Wen et al.，2017）。扬子地块的古生界富有机质泥页岩分布广泛，是我国页岩气勘探的重要目标层，在上扬子区的涪陵、长宁、昭通、威远等地页岩气田已开始工业化开采（Liu Honglin et al.，2010; Nie Haikuan et al.，2020）。但是下扬子区的海相古生界暗色泥页岩在中生代经历了强烈的构造变形和岩浆活动，不利于页岩气的保存（Pan Jiping et al.，2011; Cai Zhourong et al.，2015）。
-
前人根据地质观察和矿物温度计的计算，发现接触变质晕的宽度变化范围很大，为侵入体宽度的30%~250%（Cooper et al.，2007; Santos et al.，2009; Aghaei et al.，2015; Li Xia et al.，2016），这给深部找矿和资源潜力评价工作带来了挑战。因此，研究花岗质岩浆侵位对围岩裂隙发育和热结构的影响不仅有助于认识岩浆演化、流体迁移和区域构造的关系，而且对揭示与花岗岩相关的成矿过程及油气成藏具有重要意义。动力学数值模拟已成为研究岩浆侵位过程的有力手段，考虑到岩浆侵位方式（Galushkin，1997; Wang Dayong et al.，2011a）、岩石热参数的温度依赖性（Wang Dayong et al.，2010）、孔隙水的气化等因素（Barker et al.，1998），前人发展了不同的热传输模型，但主要关注岩浆与围岩之间的热传导，只有少数研究考虑了孔隙水的对流散热效应（Galushkin，1997; Wang Dayong et al.，2011b，2015）。对岩浆侵位导致的裂隙发育、围岩传热方式与变质晕之间的关系尚未开展模拟研究。
-
本文首先介绍了岩浆侵位导致的接触变质作用和脆性变形，然后总结了岩浆侵位的数值模拟原理并对比了不同方法的优缺点，最后使用离散元软件MatDEM建立二维岩浆侵位模型，模拟花岗质岩浆侵位导致的围岩裂隙发育过程和热结构演化，探讨影响变质晕宽度的因素以及裂隙分布与成矿的关系。
-
1 岩浆侵位过程的变质和变形
-
1.1 接触变质作用
-
接触变质作用是在侵入体周围发生的变质作用，可分为热接触变质和接触交代变质，二者共同作用形成了接触变质晕（Winter，2001）。高温不仅促进矿物的重结晶，使矿物的粒径增加，还能促进吸热化学反应形成新的矿物组合，远离侵入体时，围岩的粒径和重结晶程度迅速减小，粒径变化与温度梯度一致。与此同时，花岗质岩浆侵位过程中不断结晶分异，产生的富硅流体会渗透到围岩中形成化学梯度，热液交代围岩使其化学成分、矿物成分和结构都改变，发生接触交代变质作用。
-
受岩浆侵位的热效应、化学效应以及围岩物质组成的影响，接触变质晕可出现分带。碎屑岩变质带中常见堇青石、矽线石、红柱石等变质矿物，形成角岩带和板岩带（Brown，2013; Fu Xiaofang et al.，2014; Liu Chao et al.，2015），而碳酸盐岩变质带中常形成矽卡岩带和大理岩带（Svensen and Jamtveit，1998; Ouyang Yongpeng et al.，2018）。岩浆的类型、深度、温度、规模以及围岩的岩性和物理性质等都会影响接触变质晕的宽度，例如:泥页岩的温度敏感性高于砂岩，更易发生热变质，形成的变质晕比砂岩的变质晕更宽（Liu Chao et al.，2015）。值得注意的是，矿物温度计和有机质成熟度仅能记录围岩经历的最高温度，难以精细反映围岩的热演化历史和交代变质作用，因此用来确定接触变质晕的范围具有局限性。
-
1.2 脆性变形
-
岩浆侵位过程中会形成岩株、岩床、岩墙等不同形态的侵入体，围岩的层理、面理等软弱面是形成岩床和岩盖的有利条件，先存的高角度断裂则有利于形成岩墙。岩浆上升时对围岩的挤压和黏性拖曳会造成围岩的脆性破裂和韧性变形（Wang Tao et al.，2000; Brown，2013; Liu Chao et al.，2015）。此外，围岩中的孔隙流体受热膨胀，可在围岩中形成高孔隙流体压并使围岩发生脆性破裂，而这些破裂又会成为熔/流体迁移的通道。例如:南极洲Dronning Maud Land石英正长岩中的伟晶岩脉相互截切，延伸较远（图1a），变质晕在岩脉两侧对称分布，与未受影响围岩之间的边界清晰（图1b），变质晕内裂隙密度显著升高，变质晕的宽度对应于高渗透裂隙带的宽度，表明含水熔体导致的高孔隙流体压使围岩形成了张性裂隙，熔体沿裂隙迁移并冷却结晶（Engvik et al.，2005）。在澳大利亚西部Yilgarn克拉通Jupiter金矿，正长岩侵入绿片岩相变质玄武岩并形成不规则、放射状分布的裂隙网络，控制了含金石英-碳酸盐岩脉的分布（Duuring et al.，2000; Das et al.，2014）（图1c）。新疆阿尔泰可可托海3号脉具有典型的同心环状分带，伟晶岩脉沿着深成岩体的中等倾角节理、叶理等张性裂隙分布（Qin Kezhang et al.，2021b）。
-
Zhang Wen et al.（2017）对浙江安吉罗村花岗岩脉及其围岩的研究发现，岩脉侵入过程中的应力变化导致荷塘组页岩发育了三期微裂隙:① 岩体侵入挤压围岩，形成广泛分布的共轭剪节理（裂隙I）; ② 高温导致页岩中的黏土矿物发生脱水反应，有机质裂解成CO2及CH4，破坏了原来的纳米级有机质孔隙，高孔隙流体压和热膨胀效应导致靠近侵入体的页岩形成张裂隙（裂隙II），富硅流体沿裂隙与页岩发生交代作用，并携带碳质残留体迁移; ③ 岩体冷却后期，孔隙流体压和热收缩导致页岩发育近平行层面的收缩缝（裂隙 III）（图2a、b）。罗村花岗岩脉的侵入显著提高了页岩的有机质成熟度、裂隙密度和渗透率，从而使页岩气沿裂隙逃逸，影响范围约为侵入体宽度的37%，由于渗透率的提升可促进孔隙水对流，导致热变质晕与由裂隙密度控制的渗透率晕的范围基本一致（图2c~e）。
-
图1 岩浆侵位导致的围岩破裂和成矿作用
-
Fig.1 Fractures and mineralization of country rocks due to magma emplacement
-
（a，b）—南极洲Dronning Maud Land的伟晶岩脉及其蚀变晕（据Engvik et al.，2005）;（c）—澳大利亚Jupiter金矿（据Das et al.，2014）;（d）—新疆可可托海3号脉（据Qin Kezhang et al.，2021b）
-
(a, b) —Pegmatite vein and its alteration aureole in Dronning Maud Land, Antarctica (after Engvik et al., 2005) ; (c) —the Jupiter gold deposit in the Yilgarn Craton, western Australia (after Das et al., 2014) ; (d) —the Koktokay No.3 pegmatite (after Qin Kezhang et al., 2021b)
-
2 岩浆侵位的主要数值模拟方法
-
2.1 岩浆侵位的热传导和热对流数值模拟
-
岩浆侵入围岩涉及多个耦合的物理过程，如:热传导、裂隙生成与扩展、裂隙内熔/流体的热对流。如前所述，热接触变质作用主要由热传导控制，因此传热模型可用于定量分析侵入体导致的热效应（Jaeger，1959; Peacock，1990; Galushkin，1997; Zhang Jian et al.，1997）。Lovering（1935）根据傅里叶热传导理论推导了岩浆侵入围岩的一维热传输方程:
-
式中，T、t分别代表温度（℃）和时间（s），α为热扩散率（m2/s），x代表围岩到岩体边界的距离（m），Vx是流速（m/s），A是体积产热量（J/（s·m3）），ρ为密度（kg/m3），C是单位质量的热容（J/K）。公式（1）反映了控制围岩温度随时间变化的三个因素:右边第一项是热传导带来的烘烤效应，与围岩的热扩散率和温度的二阶导数有关; 第二项是热对流，与孔隙流体的流速和地温梯度成正比; 第三项是围岩的热源和散热，围岩的放射性产热是热源，孔隙水的蒸发、对流和矿物变质反应的吸热是主要的散热。岩石的热扩散率α（m2/s）和热导率λ（W/（m·K））遵循以下关系:
-
求解热传输方程需要给定初始条件和边界条件，包括岩体宽度、初始岩浆温度、围岩温度、热边界条件等。
-
Barker et al.（1998）研究了澳大利亚Gippsland盆地的基性岩墙侵入含碳砂岩导致的热变质晕，他们使用流体包裹体和镜质体反射率计算了含碳砂岩的古地温，将岩墙侵入导致的围岩增温模型分为四类:① 不考虑孔隙水蒸发的简单热传导模型; ② 孔隙水吸收热量并蒸发散热的复杂热传导模型; ③ 围岩快速升温，岩墙附近的高温流体上升，但未形成补给循环系统的初始对流模型; ④ 岩墙附近的高温流体上升，而远离岩墙的低温流体提供补给的对流环模型（图3）。假定D为岩墙宽度，X为围岩到岩墙边界的距离，X/D将岩墙宽度的影响归一化，从而定量研究传热方式对热变质晕宽度的影响。简单热传导模型预测的热变质晕宽度超过岩墙宽度的两倍（图3a），如果考虑孔隙水的气化吸热效应，复杂热传导模型形成的热变质晕宽度与岩墙宽度基本相当（图3b）。
-
图2 浙江安吉罗村剖面荷塘组页岩中的三期微裂隙（a，b）以及古地温（c）、裂隙密度（d）、渗透率剖面（e）（据Zhang Wen et al.，2017修改）
-
Fig.2 Three-stage fractures in shales of shales of the Hetang Formation (a, b) in Locun section, Anji, Zhejiang Province, and paleogeothermal temperature (c) , fracture density (d) and permeability profiles (e) (modified after Zhang Wen et al., 2017)
-
Raymond et al.（1988）观察苏格兰Midland Valley的拉斑质和碱性岩席侵入石炭系地层的热效应，发现拉斑质岩席的围岩压实紧密且比较干燥，主要发生热传导，产生很宽的热变质晕; 但是碱性岩席的围岩疏松且饱水，形成非常窄的热变质晕。Ma Yemu et al.（2013）计算了华南36个花岗岩侵入体-泥质围岩的二维简单热传导模型，与根据镜质体反射率估算的泥岩古地温剖面对比，发现虽然花岗岩体的出露宽度在2~30 km，但是其热影响的范围都在3 km以内。这意味着仅考虑热传导效应，难以解释侵入体的规模与接触变质晕宽度的关系。
-
由公式（1）可知，围岩中孔隙流体的流速会影响热对流传递的热量，而孔隙流体的流速受许多因素的影响，例如:围岩的热应力、孔隙率和渗透率（Hayba et al.，1997）。随着接触带持续受热，孔隙流体气化，受到浮力向上运移，导致远处的低温地层水向岩墙方向迁移，发生初始热对流，形成的热变质晕宽度显著小于岩墙宽度（图3c）; 随着加热时间的延长，宽岩墙附近可形成完整的对流环，使得变质晕的宽度增大，变质晕宽度可以达到岩墙宽度的2倍以上，并伴随温度波动（Delaney，1987; Barker et al.，1998）（图3d）。但是，Barker et al.（1998）只基于一维热传导模型计算了简单热传导和复杂热传导的热变质晕宽度，并没有对初始对流模型和对流环模型开展模拟验证。
-
大多数岩浆侵位的热传输模型假设侵入体瞬间侵入地层，围岩中的孔隙流体是否会发生对流以及在何处发生对流可以通过瑞利数（Ra）来确定。瑞利数是与热浮力驱动的对流相关的无量纲数，与流体的热力学性质和介质的渗透率密切相关（Fu et al.，2010）:
-
式中，g是重力加速度（m/s2），βf、ρf、Cf和ηf分别为流体的体积膨胀系数（K-1）、密度（kg/m3）、比热（J/（kg·K））和动力学黏度（Pa·s），λ是围岩的热导率（W/（m·K）），k为围岩的渗透率（m2），d是对流单元的厚度（m），ΔT是对流单元的温度差（K）。当流体的瑞利数低于临界值时，热量传递的主要形式是热传导; 当超过临界值时主要是热对流。Fu et al.（2010）假定围岩的渗透率为10-14~10-17 m2，获得围岩中Ra的临界值为4π2。Hayba et al.（1997）发现当围岩的渗透率小于 10-16 m2 时，热传导是围岩中主要的传热模式。由于页岩的渗透率通常低于10-16 m2（Brooks Hanson，1995; Zhang Wen et al.，2017，2018），因此前人在模拟岩浆侵入泥页岩时，经常忽略热对流的影响（Wang Dayong et al.，2007; Ma Yemu et al.，2013）。然而，基性岩床侵入页岩形成的不对称热蚀变晕表明，热对流可能在接触变质晕的发育中扮演重要角色（Simoneit et al.，1981; Qiu Longwei et al.，2000）。
-
图3 岩墙侵入后围岩的四种加热模型（a~d）
-
Fig.3 Four heating models (a~d) of host rocks after intrusion of a dike
-
D —岩墙宽度; X —围岩到岩墙的距离; X/D —到岩墙的归一化距离; Tc —侵入瞬间的接触温度; 箭头指示流体的运动方向（据Barker et al.，1998修改）
-
D —Width of a dyke; X —distance to the dyke; X/D —the normalized distance to the dyke; T c—the contact temperature at instant of intrusion; arrows indicate the motion direction of the fluids (modified after Barker et al., 1998)
-
Nu为努塞尔数，等于水热传输与仅通过热传导进行的热传输的比值，通常在0~4的范围内（Krynauw et al.，1994; Galushkin，1997）。目前主要通过将围岩的热导率提高（1+Nu）倍来估计热对流对整体热量的影响，该方法可以得到岩床周围的不对称变质晕（Galushkin，1997; Wang Dayong et al.，2011b），但无法比较不同对流模式对围岩热结构的影响。Wang Dayong and Manga（2015）使用热浮力作为对流的驱动力，通过孔隙空间计算热对流。但是，对流方式与孔隙的连通性密切相关，通过孔隙发生的渗流与裂隙主导的层流对围岩的热结构会有不一样的影响。因此，在分析岩浆侵位导致的接触变质作用时，必须考虑裂隙的性质、空间分布及热对流效应。
-
2.2 围岩裂隙发育的数值模拟
-
对于岩浆充填的裂隙扩展，驱动岩浆在裂隙中流动的应力来源包括:岩浆源区的超压，岩浆的浮力，以及垂直于侵入体边界的构造应力梯度（Rubin，1995）。连续介质法一般采用连续介质线弹性断裂力学、连续损伤力学或者结合裂隙模型的方法进行分析。例如:Das et al.（2014）使用Grady-Kipp连续损伤模型，结合光滑粒子流体力学法（Smooth particle hydrodynamics method，简称SPH）对澳大利亚西部Jupiter金矿的侵入体进行建模，通过弹性变形-流体流动-热演化耦合的模型计算围岩变形和裂隙发育，发现在构造伸展背景下模拟的裂隙分布结果与Jupiter金矿的矿脉分布特征一致。Chen Tielin et al.（2017）结合弥散裂缝模型，应用有限元法研究了浅层沉积物中的岩脉侵入过程，模拟结果与Abdelmalak et al.（2012）的物理实验结果相似。该方法将裂隙的影响等效于材料本构模型的变化，产生裂隙时只调整材料的本构矩阵，从而避免连续介质网格划分的不足，但是无法考虑裂缝的实际分布及局部应力状态。
-
非连续介质法将介质看成离散的几何单元，利用离散粒子之间的相互作用来定义系统的聚合行为。在离散元法中，岩土体模型由一系列遵循牛顿运动定律的颗粒堆积组成，对每个离散的几何体颗粒都赋予材料特性（如密度、弹性模量等）和颗粒间特性（摩擦、剪切和拉伸强度等），基于牛顿第二定律计算颗粒的加速度、速度和位移。在基本的线弹性接触模型中，颗粒间通过弹簧力相互作用（图4），两颗粒间的法向力（Fn）由下式给出（Liu Chun et al.，2017）:
-
式中，Fn为法向力（kN），Kn为法向刚度（kN/m），Xn为法向相对位移（m），Xb为断裂位移（m）。初始时，颗粒与其相邻颗粒相互连接，受拉力或压力作用（公式4a），当颗粒之间相对位移超过断裂位移时，弹簧断裂，颗粒间拉力消失（公式4c），仅存在压力作用（公式4b）。
-
弹簧间通过切向弹簧来模拟剪切变形和切向力，剪切弹簧的切向力Fs为:
-
式中，Ks为切向刚度（kN/m），Xs为切向相对位移（m）。对于完整的颗粒连接，其最大切向力Fsmax由摩尔-库仑岩石破裂准则确定:
-
式中，Fs0为颗粒内部的抗剪力（kN），μp为摩擦系数，当颗粒受切向外力超过最大切向力Fsmax时，颗粒间的连接断裂。
-
图4 离散颗粒作用示意图（据Zhang Hongyong et al.，2021）
-
Fig.4 Schematic diagram of the action of discrete particle (after Zhang Hongyong et al., 2021)
-
（a）—离散颗粒堆积模型;（b）—颗粒间法向弹簧力;（c）—颗粒间切向弹簧力
-
(a) —Discrete particle accumulation model; (b) —normal spring force between particles; (c) —tangential spring force between particles
-
与连续介质法相比，离散元法最大的优势是可以模拟真实的不连续面，因此适合模拟破碎过程、裂隙和断层演化等。Zhao Chongbin et al.（2008）利用离散元颗粒流程序（Particle flow code，简称PFC）模拟岩盖侵位导致上覆围岩形成穹隆的过程，他们采用可压缩性相对较小的流体粒子来模拟岩浆，用常规固体粒子模拟围岩，当粒子间的接触力超过粒子间的拉伸黏结强度，形成张裂，当剪切接触力超过相应的剪切黏结强度，则产生剪裂，裂隙进一步扩展连通形成断裂。但是，离散元法在地质上的应用还面临许多挑战，主要包括:建模困难，材料宏观性质和颗粒微观特性之间的转换关系尚不明确; 计算量非常大，在大尺度模拟中观察微小变形很困难; 岩土体由颗粒、孔隙和孔隙中的流体组成，流-固耦合问题是离散元模拟中的难点（Gray et al.，2014; Liu Chun et al.，2020）。此外，离散元法模拟的参数与自然界的地质参数难以直接对比，因此前人借鉴构造物理模拟实验的相似性原理，对离散元法模型进行放缩以模拟大尺度、长期演化的地质过程（Feng Yuntian et al.，2014; Zhou Zongqing et al.，2019; Li Changsheng et al.，2021）。
-
2.3 孔隙密度流法模拟流-固耦合
-
针对非连续介质的流-固耦合问题，南京大学刘春教授及其团队基于矩阵离散元软件（Matrix Discrete Element Method，简称MatDEM），采用孔隙密度流法实现离散元中的流-固耦合数值模拟，成功用于隧道注浆、海床沉积物孔压的累积模拟（Wang Yue et al.，2021; Zhang Hongyong et al.，2021）。该方法建模的基本步骤是:首先堆积颗粒（图5a），然后根据颗粒间的接触关系，通过连接相邻颗粒的中心点，在堆积模型中剖分出一系列相互连通的流体域（图5b），最后得到孔隙流体网络-固体颗粒骨架模型（图5c）。假定孔隙水饱和，孔隙流体压力Pf由孔隙流体的密度ρf和温度T确定（Liu Chun et al.，2020）:
-
图5 单元堆积和流体域示意图（据Zhang Hongyong et al.，2021）
-
Fig.5 Schematic diagram of unit accumulation and fluid domain (after Zhang Hongyong et al., 2021)
-
（a）—离散元堆积模型;（b）—孔隙流体域;（c）—颗粒单元和孔隙系统
-
(a) —Discrete element stacking model; (b) —pore fluid domain; (c) —particle unit and pore system
-
当相邻孔隙间存在压力差时，流体将通过二者间的孔喉通道渗流（图6），其渗流量可用类似达西定律的方法来计算（Wang Yue et al.，2021）。模型中单位时间内通过孔喉的流量q定义为:
-
式中，dP为孔喉间压力差（Pa），定义为孔喉的渗透系数（m/s）。K为孔喉的渗透系数因子（s-1），类似于岩石的宏观渗透系数，与颗粒的形状、大小、排列方式以及流体的密度和动力学黏度有关，通过修改K可以改变围岩的渗透系数，但岩石的宏观渗透系数与微观的孔喉渗透系数因子K之间的转换关系尚未进行标定。B为孔喉通道面积（m2），孔喉长度l（m）定义为相邻颗粒中较小的颗粒半径:
-
式中，R1、R2是构成孔喉的两颗粒半径（m）。对于二维问题，单位宽度孔喉通道面积可用孔喉的直径dw表示。由于二维颗粒堆积通常会封闭孔喉，使得孔喉直径为0，无法发生渗流，因此需要给每个单元定义较小的渗流作用半径 Rw，默认取单元半径的0.975倍，从而保证一定的孔喉直径，使流体能够通过孔喉运移。
-
孔喉的直径dw定义为:
-
式中，L为颗粒中心距离（m），Rw1和Rw2分别为两个单元的渗流作用半径（m）。图6中的矩形区域定义了孔隙P1和P2间的渗流通道。以上计算中，流体在孔隙压力差的驱动下运移，发生渗流。模拟时流体压力作用于固体颗粒引起颗粒运动，而固体颗粒的位移变化又导致相邻孔隙的体积、密度和压力变化，引起渗流，从而实现流-固耦合作用的模拟。
-
图6 颗粒微观渗流示意图（据Wang Yue et al.，2021）
-
Fig.6 Schematic diagram of particle micro seepage (after Wang Yue et al., 2021)
-
3 岩浆侵位的离散元建模
-
为研究花岗质岩浆侵位过程中的流-固耦合和接触变质晕范围，本文利用MatDEM建立了双层围岩和均质围岩两种模型。双层围岩模型用于验证MatDEM软件模拟裂隙发育的可靠性，探索岩浆侵位过程中裂隙的发育过程; 均质围岩模型用于计算岩浆冷却过程中，不同传热方式对变质晕范围的影响（图7）。本文通过对MatDEM软件的进一步开发，模拟岩浆侵位过程中的流-固-热-力耦合。
-
对于双层围岩模型，参考Das et al.（2014）的模型，首先建立长4400 m、深1200 m的二维模型箱，向内随机添加平均半径为11 m的颗粒，总颗粒数为12049个，经过重力沉积和压实得到颗粒堆积模型（图7a）。对于双层围岩模型，上部200 m为自由空气层，纵坐标Z的原点为空气层的上表面，因此，距离地表的深度等于当前点的Z坐标减去200 m。堆积完成后，删去模型箱两侧各200 m范围内的颗粒，得到长4000 m、深1200 m的地质模型。
-
根据Das et al.（2014）的材料力学性质（表1），基于宏微观转换公式（Liu Chun et al.，2017）和MatDEM自动训练材料模块，得到符合宏观性质的离散元微观力学参数（表2）。根据构造物理模拟实验中的相似性理论，Li Changsheng et al.（2021）对比了用离散元法建立的不同尺度的模型，发现不同尺度下的模型变形是等效的。因此，本文根据等比缩放原则，对模型的离散元微观力学参数进行缩放，重力加速度降为原来的10-7，材料为各向同性，密度增大为原来的108倍，对应的模拟时间增大1000倍，由于孔喉渗透系数因子K与流体黏度相关，因此黏度的缩放包含在岩石宏观渗透系数的缩放中。
-
如图7a所示，双层围岩模型的上层为绿岩，下层为花岗岩。将绿岩底面与花岗岩接触的颗粒间的抗拉强度和抗压强度降为原来的1/5，模拟绿岩和花岗岩的岩性界面。删除花岗岩中椭圆形区域的颗粒，建立流体域，作为岩浆房（图7a）。在岩浆注入点保持50 MPa侵位压力，模拟岩浆持续注入。地面设置恒定大气压0.1 MPa，地表温度为25℃，岩浆的温度恒定为727℃（~1000 K），地温梯度30 km/℃，不考虑围岩的放射性生热、岩浆结晶释放的潜热以及流体气化吸收的热量。根据含水量5%的花岗质岩浆在~1000 K时密度随压力的变化（Malfait et al.，2014），拟合得到岩浆的密度-压力曲线。将系统中的流体状态方程参数修改为岩浆状态方程的参数，从而得到岩浆在孔隙中迁移时压力和密度的变化关系。为减少模拟时间，本研究使用水的黏度来代替花岗质岩浆的黏度，且不考虑温度对黏度的影响，即模拟过程中体系内熔体的黏度不变。围岩的宏观渗透系数与微观上的孔喉渗透系数有关，通过修改公式（8）中的孔喉渗透系数因子来设定围岩的渗透系数。
-
图7 岩浆侵位的初始模型设置
-
Fig.7 Model setting for magma intrusion
-
（a）—绿岩和花岗岩组成的双层模型;（b）—岩墙模型;（c）—穹隆模型
-
(a) —A double-layer model composed of greenstone and granite; (b) —a dike model; (c) —a dome model
-
模型箱的下边界在垂直方向固定，水平方向可以平移; 左右边界在垂直方向和水平方向上都可以移动。左右边界的水平运动速率分别为0.01 m/s，模拟构造伸展的地质背景。模型计算的时间步为5×10-4 s，每个周期迭代计算1000次，加载周期0.5 s，循环加载20次，模拟时间为0.5 s×20次=10 s。根据相似性原理进行放大，则对应自然界中岩浆侵位时间为10 s×1000=104 s，相当于40 km×10 km的二维地质体中岩浆侵位2.7 h的演化结果。
-
在模型中，热量被认为储存在流体中。热传导发生在两个温度不同的相邻流体区域之间，热量从高温流体域自发地转移到低温流体域。热传导的方向与温度梯度的方向一致，热传导效率随着温度梯度的增加而增加。热对流则是通过孔隙和微裂隙之间流动的流体携带热量，导致温度变化。热对流的传播方向与压力梯度一致，热对流效应与相邻孔隙之间的压力差有关，压力差越大，流速越快。在双层模型中，随着岩浆持续注入，颗粒间的接触力变化可在围岩中形成张裂隙或剪裂隙，导致孔喉的渗透系数增加。如果裂隙与岩浆房连通，岩浆就会沿裂隙贯入，裂隙中熔体的温度与岩浆温度相同。因此，一方面岩浆与围岩之间不断发生热传导，另一方面岩浆房的压力促使熔体在围岩的孔隙空间及裂隙中流动，通过对流传递热量。
-
为对比不同热传输机制对变质晕的影响，本文使用均质围岩模型计算岩浆侵入导致的围岩增温过程，即岩浆的热烘烤效应。模型假设:① 岩浆瞬间侵位，侵位后不流动; ② 围岩为均质体; ③ 不考虑岩浆挥发分释放引起的能量及质量损失; ④ 围岩热导率、流体热膨胀系数为常数，不随温度变化; ⑤ 忽略孔隙水蒸发散热; ⑥ 不考虑围岩和岩浆的放射性生热。在MatDEM中建立长1000 m，深500 m模型箱，颗粒平均半径4 m，总颗粒数8656个。赋予绿岩的材料力学参数，岩体的力学性质和热传导率为各向同性，模型中建立矩型的流体域作为岩墙（图7b），不设置构造应力，岩浆侵位温度为727℃，温度边界条件与双层模型相同。加热时围岩孔隙中的水受热产生热增压，即孔隙压随孔隙水的温度增加而增加，热增压系数由水的等压体积热膨胀系数和等温体积压缩系数的比值决定，热增压驱动流体发生对流。
-
4 模拟结果
-
4.1 双层围岩模型的裂隙发育
-
在双层模型中，高温岩浆持续侵入，诱发裂隙，围岩中不仅存在热传导，还通过裂隙发生热对流。图8为构造伸展背景下迭代20000次时颗粒的位移结果。在颗粒的垂向位移场中，受岩浆浮力作用，岩浆房整体向上推挤导致上方颗粒向上运动，如图8a中的红色区。在岩浆房斜上方约60°的方位，颗粒间的相对位移超过断裂位移，弹簧断裂，颗粒间拉力消失，形成裂隙，表现为颗粒垂向位移场的快速变化（图8a）。在颗粒的水平位移场中，岩浆房侧向推挤围岩，在岩浆房两侧形成位移方向相反、近对称分布的最大水平位移区，但是并没有形成水平方向展布的裂隙（图8b）。结合两个方向的位移场，可以看到在岩浆侵位的过程中，岩浆房上部发育了一组近似对称分布的主干张性断裂，两条主干断裂的锐夹角以内的围岩向上运动，导致垂向位移场中地表轻微抬升。
-
图9a~d为双层模型在构造伸展背景下迭代2500次、5000次、7500次、20000次时围岩中裂隙的发育情况。侵位初期受岩浆的挤压作用，在岩浆房上方与围岩接触区发育裂隙，并不断向远处延伸（图9a、b），大体呈径向展布，同时岩浆的挤压使得岩浆房上方绿岩和花岗岩的接触界面形成高应力集中，导致裂隙沿界面扩展。迭代7500次时，大量的径向裂隙已扩展至地表（图9c），绿岩和花岗岩的界面裂隙也不断向模型两侧延伸。迭代20000次时，随着岩浆不断侵入，岩浆房与围岩接触区的孔隙流体压不断增加，使得此区域内发育大量张性裂隙，裂隙密度急剧增大，模型最终形成了三个高裂隙密度区:绿岩和花岗岩的岩性界面、岩浆房上方围岩、岩浆房与围岩的接触边界（图9d）。
-
本模型中颗粒的力学性质根据相似性原理缩放可以与真实世界对比，而摩尔-库仑岩石破裂准则不依赖于应变率，模拟的时间与变形事件的持续时间没有直接关系。由于使用水的黏度代替岩浆黏度，模拟的裂隙形成和扩展过程比天然岩浆侵位时的裂隙扩展更快。本模型的裂隙发育过程可以较好地与Das et al.（2014）使用光滑粒子流体力学法计算的142 ka裂隙发育模型（图9e~h）对比，验证了使用MatDEM孔隙密度流法模拟岩浆长期侵位过程的可行性。但是，两个模型给出的裂隙发育过程有所不同，将在下文具体讨论。
-
图8 岩浆侵位后双层模型中颗粒的位移场
-
Fig.8 Displacement distribution of particles due to magma intrusion in a double-layer model
-
（a）—垂直方向的位移，红色为向上运动，蓝色为向下运动;（b）—水平方向的位移，红色为向右运动，蓝色为向左运动
-
(a) —Displacement in vertical direction, red color for upward motion and blue color for downward motion; (b) —displacement in horizontal direction, red and blue colors infer to motion towards right and left, respectively
-
图9 构造伸展下围岩中的裂隙发育模式对比
-
Fig.9 Comparison of fracture development patterns in host rocks under tectonic extension
-
（a~d）—本文MatDEM双层模型的模拟结果;（e~h）—Das et al.（2014）的模拟结果
-
(a~d) —Simulation results of a double-layer model using MatDEM in this study; (e~h) —simulation results of Das et al. (2014)
-
图10为双层模型中的温度场变化，其中红色代表岩浆熔体，温度为727℃，从岩浆房到围岩形成温度梯度（图10a）。熔体在浮力和压实作用下的运移方式可归纳为孔隙流和通道流两种，孔隙流发生在颗粒尺度上，可以遍布整个模型。结合裂隙发育过程（图9）可知，岩浆侵位时，靠近侵入体的围岩中产生大量微裂隙，当微裂隙连通，形成断裂并与岩浆房连通（图10b），岩浆以通道流的形式在贯通的裂隙中快速迁移，冷却后形成岩脉（图10c）。左侧裂隙扩展更快并抵达地表，因此岩浆沿着贯通的张裂一直侵位到地表附近，在水平拉张作用力下，岩浆由斜向上侵位逐渐转为近垂直侵位（图10d）。随着岩浆侵位时间延长，岩浆房周围变质晕的宽度不断增加，而岩脉的形成进一步增加了变质晕的影响范围。值得注意的是，虽然在绿岩和花岗岩的界面也形成了大量微裂隙，但是岩浆并没有沿着岩性界面侵入形成岩席，将在下文具体讨论。
-
4.2 均质围岩模型的热变质晕
-
本文使用均质围岩模型定量研究热对流和侵入体形状对热变质晕的影响（表3）。设置流体域间的导热系数不变，改变围岩的孔喉渗透系数因子和传热方式得到三组对比实验1~3。此外，由于花岗质岩浆黏度较大，往往形成穹隆状侵入体，因此修改岩墙侵入体为椭球状的穹隆侵入体（图7c），保持流体域间的导热系数不变，设置围岩不同渗透系数和不同传热方式的对比实验4~6。为探究裂隙发育对围岩热传导模式的影响，向岩浆房施加恒定温度和孔隙流体压，使其持续侵位并在围岩中形成裂隙，将初始侵位阶段的含裂隙围岩模型用于计算热传输过程（表3中的实验7）。
-
假定侵入后的岩墙原位冷却，实验1仅考虑热传导，实验2同时考虑热传导和热对流，围岩的渗透率不变，对围岩到岩墙的距离进行归一化。从同一时刻的围岩温度剖面来看，热对流的参与导致岩墙附近的温度更高，接触带围岩的温度梯度更陡（图11）。在336 s，只发生热传导的围岩已经基本达到温度平衡（图11a），而同时发生热传导和热对流的情况下，在侵入体边界的围岩温度仍高于远处的围岩温度（图11b），表明热对流提高了侵入体附近围岩的温度，减缓了接触带冷却的速率。
-
图10 岩浆侵位阶段模型的温度场演化（a~d）
-
Fig.10 Evolution of temperature field (a~d) during magma intrusion in a double-layer model
-
假定200℃为热变质晕的温度下限，图12为相同渗透率的围岩在岩墙侵入后，在热传导（实验1）和热传导加热对流（实验2）两种模式下，变质作用范围最宽时的温度场分布，即热变质晕范围，以便与使用地质温度计获得的变质晕对比，图中超过200℃的围岩都显示为红色。以岩墙边界的中点为起点，截取围岩中水平方向和垂直方向的温度剖面，获得围岩到侵入体不同距离处的峰值温度（Tpeak）曲线，峰值温度达到200℃的位置即变质晕的范围。对于直立岩墙，水平方向的变质晕宽度（剖面1）显著大于垂直方向（剖面2）。与仅通过热传导升温的围岩（图12a、b）对比，当热对流参与热量传输时，侵入体附近围岩的峰值温度略微升高，但归一化的变质晕宽度X/D减小（图12d），意味着热对流使岩浆的热烘烤效应更集中，这与不同传热机制形成的某一时刻的围岩温度剖面模拟结果（图11）一致。
-
由公式（8）可知，孔喉的流量与岩石的渗透系数成正比，并控制了岩石中热对流的效率。为研究渗透系数对接触变质晕发育的影响，实验3假定均质围岩的孔喉渗透系数因子为10-5，计算热对流参与情况下的变质晕达到最宽时的温度场分布。与实验2围岩的孔喉渗透系数因子为10-6对比（图12d），渗透系数升高导致接触变质晕的范围显著变窄，意味着围岩的渗透系数越大，热对流效应越强，热增压效应导致远处围岩中的低温流体快速向岩墙方向流动，使靠近岩墙的围岩快速冷却（图13）。
-
图11 直立岩墙侵入形成的均质围岩温度剖面
-
Fig.11 Temperature profiles of homogeneous host rocks due to intrusion of a vertical dike
-
（a）—仅考虑热传导;（b）—热传导和热对流; D —岩墙宽度; X —围岩到岩墙的距离; X/D —到岩墙的归一化距离
-
(a) —Only heat conduction; (b) —heat conduction and convection; D —the width of a dyke; X —the distance to the dyke; X/D —the normalized distance to the dyke
-
图12 直立岩墙侵入均质围岩形成的接触变质晕
-
Fig.12 Contact aureole in homogeneous host rocks due to intrusion of a vertical dike
-
（a、b）—仅考虑热传导;（c、d）—热传导和热对流; 白色虚线为温度剖面位置
-
(a, b) —Only heat conduction; (c, d) —heat conduction and convection; white dashed lines indicate the positions of two temperature profiles
-
对于原位冷却的穹隆状岩浆房，均质围岩模型的实验4仅考虑热传导，围岩的渗透系数因子为10-6，实验5和实验6同时考虑热传导和热对流，围岩的渗透系数因子分别为10-6和10-7。截取围岩中水平方向和垂直方向的温度剖面，仍使用岩墙宽度对围岩到岩墙的距离进行归一化。热传导形成的变质晕范围最宽（图14a、b），热对流降低了变质晕的宽度（图14c、d），与岩墙的模拟结果一致。实验6的结果表明:当围岩的渗透系数极低时，对流难以发生，可以忽略热对流的影响，围岩的传热方式以热传导为主（图14e、f）。对比岩墙冷却和穹隆冷却模型，可知侵入体的几何形态控制了接触变质晕的空间形态。
-
4.3 岩浆侵位初期的裂隙发育和热传导
-
当岩浆侵入到均质围岩中，岩浆将挤压围岩。如图15a所示，对于穹隆状花岗质岩浆房，侵位前围岩的初始平均孔喉渗透系数为4.0663×10-6，在迭代500步（约0.5 s）时降低到2.7705×10-6，这表明岩浆开始侵位时，岩浆的挤压力使围岩中的软孔隙迅速关闭，减低了围岩的孔喉通道面积（公式8中的B），从而导致孔喉渗透系数降低。随着岩浆的持续侵位，围岩的平均孔喉渗透系数逐渐缓慢增大，随时间的波动反映了微裂隙的产生（图15a）。但是直到22 s，围岩的平均孔喉渗透系数只达到2.95×10-6，仍然小于侵位前围岩的初始平均孔喉渗透系数。
-
图13 热传导和热对流共同作用下直立岩墙侵入均质围岩形成的接触变质晕
-
Fig.13 Contact aureole in homogeneous host rocks due to heat conduction and convection after intrusion of a vertical dike
-
图14 花岗质穹隆侵入均质围岩形成的接触变质晕
-
Fig.14 Contact aureole in homogeneous host rocks due to intrusion of a granitic dome
-
（a、b）—仅热传导，围岩渗透系数因子K=10-6;（c、d）—热传导和热对流，K=10-6;（e、f）—热传导和热对流，K=10-7
-
(a, b) —Only conduction, permeability coefficient factor of host rocks K=10-6; (c, d) —conduction and convection, K=10-6; (e, f) —conduction and convection, K=10-7
-
在侵位6 s时，在岩浆房上部的围岩中出现较多裂隙（图15b），采用此时的裂隙模型来计算围岩传热过程（表3中的实验7），与不含裂隙的围岩传热实验（表3中的实验5）对比，二者都既有热传导，又有热对流，实验7和实验5的平均孔喉渗透系数分别为2.84×10-6和4.0663×10-6。穹隆状花岗质岩浆房侵位形成的含裂隙围岩的变质晕范围（实验7，图16b）与不含裂隙围岩的变质晕范围（实验5，图14d）几乎一致。结合岩浆侵位初期平均孔喉渗透系数的变化（图15a），可知岩浆侵位初期虽然在围岩中产生了裂隙，但是在岩浆的挤压应力下，这些早期剪裂并不能成为流体迁移的有效通道，对围岩的变质晕范围影响很小。
-
5 讨论
-
5.1 裂隙发育和热传输模式
-
虽然本文的双层围岩模型获得的裂隙分布可以与Das et al.（2014）使用光滑粒子流体力学法模拟的142 ka的裂隙分布对比（图9），但是二者反映的裂隙发育机制和过程并不相同。导致差异的原因是光滑粒子流体力学法本质上仍属于连续介质方法，该方法虽然不需要划分网格，将连续的物质通过相互作用的粒子来描述，但孔隙和微裂隙相对于颗粒来说是非连续的，因此该方法无法模拟孔隙流体压的作用。在Das et al.（2014）的模型中，岩浆侵位的挤压应力使围岩发生强烈的径向压缩，导致岩石局部压碎破坏（图9e、f）。随着岩浆持续侵入，高应力区的影响范围不断增加直至顶部自由面，裂隙扩展并形成大量径向断裂（图9g、h）。而本文的双层围岩模型模拟结果较好地再现了安吉罗村剖面页岩的裂隙发育过程:岩浆侵位初期围岩中发育大量呈径向展布的剪裂隙（图9a、b），围岩中的热传导和孔隙渗流占主导（图10a、b），与安吉罗村页岩中广泛发育的第一期共轭剪裂相对应（图2a）（Zhang Wen et al.，2017）。随后在持续的孔隙流体压下，径向剪裂隙连通形成主干张性断裂并向上扩展至地表，成为熔体迁移的重要通道（图9c），通道流和局部的热对流控制了围岩的热传输（图10c）; 岩浆侵位后期，在侵入体与围岩的接触区的裂隙密度急剧增大（图9d），使热对流加强（图10d），对应于罗村页岩变质晕范围内的第二期张裂隙（图2a）。因此，本文的模拟结果解释了罗村剖面热变质晕宽度与由裂隙密度控制的渗透率晕的范围基本一致的原因（图2c、d），为评价岩浆活动对油气藏保存的影响提供了新的手段。
-
图15 花岗质穹隆初始侵位过程中围岩的平均孔喉渗透系数变化（a）及侵位6 s时的裂隙发育（b）
-
Fig.15 Change of average permeability coefficient of pore throats during initial intrusion of a grantic dome (a) and fracture development at 6 seconds (b)
-
图16 花岗质穹隆初始侵位过程中发育裂隙的围岩形成的接触变质晕
-
Fig.16 Contact aurole in fracture-bearing homogenous host rocks during initial intrusion of a grantic dome
-
值得注意的是，岩浆侵位除了需要空间，还受控于岩浆与围岩的密度差。岩浆侵位时不断冷却，与围岩的密度差逐渐减小，当密度差不足以提供岩浆继续上升的动力，而底部还有岩浆不断注入时，岩浆才会在高孔隙流体压作用下由向上侵位变为顺层横向扩展，形成岩席（Rubin，1995）。Das et al.（2014）模型中，随着裂隙的发育，岩浆沿着绿岩和花岗岩的岩性界面侵入形成岩席（图9h）。而在本文的双层模型中虽然在岩性界面的裂隙密度增加（图9d），但岩浆并没有沿着界面侵入（图10），可能是由于主干张性断裂与地表连通，提供了一个持续的低压区，从而使注入的岩浆持续上侵。本文的双层模型在5000 steps时，在侵入体上方的围岩中形成对称分布的弥散性裂隙，但是演化到7500 steps时，岩浆房上方左侧的主干断裂更为发育。我们推测这是由于控制岩石破裂的方程是非线性的，非常小的数值扰动可能被逐步放大，使左侧的主干断裂在孔隙压作用下不断增强，最终在侵入体上方形成非对称的断裂分布和岩脉分布。
-
前人研究表明，花岗伟晶岩型矿床常沿张性断裂分布，例如:阿尔泰可可托海3号脉（Qin Kezhang et al.，2021b），甲基卡锂矿化伟晶岩脉（Fu Xiaofang et al.，2014; Xu Zhiqin et al.，2018）。而热液型矿床常在花岗质岩体的接触变质晕中以脉状、网状和细脉浸染状产出，表明成矿流体沿着多期裂隙迁移、充填和沉淀（Chen Yanjing et al.，2007; Ni Pei et al.，2020）。本文的模拟结果表明:主干张性断裂提供了岩浆侵位的有利通道和伟晶岩型矿床的成矿空间，而岩浆侵位早期在岩体上方形成的弥散性剪裂在孔隙流体压的作用下被改造成张裂，有利于热液型矿床的形成（图9a~c）。江西朱溪矽卡岩型钨矿主要产于花岗岩体与碳酸盐岩接触带、滑脱构造面之上的碳酸盐岩地层中（Ouyang Yongpeng et al.，2018）。双层围岩模型的模拟结果表明，当围岩渗透率为低到中等时，在岩浆侵位后期会围绕岩浆房形成高裂隙密度区，并沿岩性界面（薄弱面）发育大量裂隙（图9d）。这些张性裂隙促进了花岗质岩浆与碳酸盐岩之间的熔-流体活动和交代作用，有助于形成矽卡岩型矿床。
-
根据数值模拟研究，岩浆侵位过程中的岩浆黏度、岩浆注入速率、围岩的杨氏模量和厚度、岩性界面的宽度和强度是控制围岩应力场的主要参数，会影响裂隙的发育过程（Giudicepietro et al.，2016; Feng Jianwei et al.，2019; Cui Zhuang et al.，2022）。高岩浆黏度、高强度围岩会产生较大的应力值，而低岩浆黏度、低强度围岩导致较低的应力，有利于岩席扩展。本文使用水的黏度代替花岗质岩浆黏度，而且没有考虑温度对黏度的影响。流体的黏度随温度升高而降低，黏度会影响瑞利数和孔喉渗透系数因子。根据公式（3），黏度越大，瑞利数越低，越难以发生热对流。根据公式（8），黏度越大，岩石的渗透系数越低，孔喉的流量越小。因此，该假设除了提高岩浆流动速率从而缩短模拟时间外，也导致围岩的应力场偏低，所以双层围岩模型在岩浆注入点采用了50 MPa的侵位压力来提高围岩的应力场。此外，本文只模拟了伸展构造背景下围岩中的裂隙发育，而斑岩型矿床往往产于汇聚板块边界，受到俯冲板片脱水引发的岛弧岩浆作用、区域性断裂等多种因素控制（Hou Zengqian et al.，2003; Chen Yanjing，2010; Zhao Junxing et al.，2012; Richards，2015），我国“硬岩型”大型锂矿也产于松潘-甘孜造山带中（Xu Zhiqin et al.，2018）。根据Das et al.（2014）的模拟结果，在拉张应力场下侵入体周围易形成放射状裂隙，而在挤压应力场下侵入体周围形成近水平裂隙。因此，后续研究需要对上述参数开展系统的模拟实验，从而深入理解不同构造背景下的花岗岩与成矿作用。
-
5.2 变质晕的宽度
-
前人的数值模拟研究表明:围岩的热导率和渗透率、对流单元的发展、侵入体的规模、岩浆温度等是影响接触变质晕宽度的重要因素（Wang Dayong et al.，2010; Wang Min et al.，2010; Wang Man et al.，2012; Wang Dayong and Manga，2015）。一般认为围岩的渗透系数较小时，围岩以热传导为主（Hayba and Ingebritsen，1997），渗透系数升高可引起初始对流，减小变质晕的范围，但当渗透系数进一步升高，远处的低温流体参与对流，出现对流环，将使变质晕范围扩大（图3; Barker et al.，1998）。
-
本文的模拟结果表明:渗透系数低时围岩以热传导为主，渗透系数升高时围岩中发生对流，变质晕的宽度变窄，符合Barker et al.（1998）根据地质观察提出的初始对流模型。通过对比模型的温度场和应力场，发现岩浆侵位初期接触带的孔隙水受热产生热增压，驱动孔隙水流向远处; 接触带的孔隙水在迅速达到峰值温度后，随着温度降低，压力逐渐降低，而远处围岩中相对低温的孔隙水不断加热，在热增压的作用下流向侵入体，导致接触变质晕变窄。这一现象符合岩浆侵入含水围岩导致的地下水流动规律（Delaney，1982）。值得注意的是，在岩墙模型中，实验3比实验2的围岩渗透系数增加了一个数量级，但并没有出现Barker et al.（1998）预测的对流环以及变质晕范围扩大的趋势，这可能是因为本文的模拟时间较短，对流环模型还有待进一步验证。
-
本文对花岗质岩墙和穹隆的模拟结果表明，虽然接触变质晕的形态受侵入体形状控制，但变质晕宽度并不一致，岩墙在水平方向的变质晕宽度大于垂直方向（图12a、b），说明岩墙在垂直方向上的热扩散比水平方向更快。这与Zhang Jian and Shi Yaolin（1997）使用二维有限元法模拟玄武质岩墙侵入沉积盆地的结果一致，虽然他们忽略了对流的影响。微裂隙会显著增强页岩的渗透率及其各向异性（Zhang Wen et al.，2018）。根据公式（8），渗透率越大，该方向孔喉的流量越大，热对流效应越显著。因此，估算矿体的成矿潜力时，必须考虑围岩的裂隙密度、产状、热传输方式以及变质晕的空间方位等因素。裂隙发育对围岩渗透率各向异性及热对流的影响将成为今后的研究重点。
-
6 结论
-
本文首次使用离散元软件MatDEM孔隙密度流法模拟花岗质岩浆侵位和冷却过程中围岩的裂隙发育和热演化，结果表明该方法可以简洁地处理岩浆侵位过程中的流-固-热-力的耦合，主要结论如下:
-
（1）岩浆侵位初期围岩受到岩浆的挤压，产生广泛分布的、呈径向展布的剪裂隙，围岩以热传导和孔隙渗流为主导; 在持续的孔隙流体压下，径向裂隙连通形成主干张性断裂并向上扩展，成为熔体迁移的通道，通道流和局部的热对流控制了围岩的热传输; 在岩浆侵位后期，岩浆房附近围岩中的孔隙流体压增大，在侵入体与围岩的接触区形成大量张裂，加强了接触带的热对流。
-
（2）热传导和热对流的共同作用会减小接触变质晕的宽度，符合前人提出的初始对流模型。当围岩的渗透系数极低时，可以忽略热对流，而渗透系数越大，热对流效应越强，变质晕越窄。接触变质晕的几何形态受侵入体的形态控制，但是变质晕宽度在侵入体的不同位置有显著差异。
-
（3）在岩浆侵位初期产生的剪裂隙虽然提高了围岩的渗透系数，但此时岩浆对孔喉的挤压作用占据主导，导致围岩整体的渗透率变化不大，这些早期剪裂并不能成为流体迁移的有效通道。
-
（4）围岩裂隙发育的模拟结果可以与地质观察对比，解释了主干张性断裂和裂隙分布与花岗伟晶岩型、热液型和矽卡岩型矿床成矿过程的相关性。
-
致谢:感谢张培震院士和刘俊来教授约稿，建模过程中得到唐春安教授、朱遥和刘辉的帮助，感谢陈衍景教授和李忠海教授提出的宝贵意见。
-
参考文献
-
Abdelmalak M M, Mourgues R, Galland O, Bureau D. 2012. Fracture mode analysis and related surface deformation during dyke intrusion: results from 2D experimental modelling. Earth and Planetary Science Letters, 359-360: 93~105.
-
Aghaei H, Gurba L W, Ward C R, Hall M, Mahmud S A. 2015. Effects of igneous intrusions on thermal maturity of carbonaceous fluvial sediments: a case study of the Early Cretaceous Strzelecki Group in west Gippsland, Victoria, Australia. International Journal of Coal Geology, 152: 68~77.
-
Agirrezabala L M, Permanyer A, Suárez-Ruiz I, Dorronsoro C. 2014. Contact metamorphism of organic-rich mudstones and carbon release around a magmatic sill in the Basque-Cantabrian basin, western Pyrenees. Organic Geochemistry, 69: 26~35.
-
Barker C E, Bone Y, Lewan M D. 1998. Fluid inclusion and vitrinite-reflectance geothermometry compared to heat-flow models of maximum paleotemperature next to dikes, western onshore Gippsland basin, Australia. International Journal of Coal Geology, 37(1-2): 73~111.
-
Brooks Hanson R. 1995. The hydrodynamics of contact metamorphism. Geological Society of America Bulletin, 107(5): 595~611.
-
Brown M. 2013. Granite: from genesis to emplacement. Geological Society of America Bulletin, 125(7-8): 1079~1113.
-
Brown M, Solar G S. 1999. The mechanism of ascent and emplacement of granite magma during transpression: asyntectonic granite paradigm. Tectonophysics, 312(1): 1~33.
-
Cai Zhourong, Xia Bin, Huang Qiangtai, Wan Zhifeng. 2015. Comparative study of the tectonic setting on the formation and preservation of paleozoic shale gas between the Upper Yangtze and the Lower Yangtze platforms. Natural Gas Geoscience, 26(8): 1446~1454 (in Chinese with English abstract).
-
Camus F, Dilles J H. 2001. A special issue devoted to porphyry copper deposits of Northern Chile. Economic Geology, 96(2): 233~237.
-
Chen Jun, Wang Rucheng, Zhu Jinchu, Lu Jianjun, Ma Dongsheng. 2013. Multiple-aged granitoids and related tungsten-tin mineralization in the Nanling range, South China. Science China: Earth Science, 56: 2045~2055 (in Chinese).
-
Chen Tielin, Cheng Shaozhen, Fang Qian, Zhou Cheng. 2017. Numerical modeling of shallow magma intrusions with finite element method. Journal of Volcanology and Geothermal Research, 333-334: 53~65.
-
Chen Yanjing. 2010. Indosinian tectonic setting, magmatism and metallogenesis in Qinling orogen, central China. Geology of China, 37(4): 854~865 (in Chinese with English abstract).
-
Chen Yanjing, Ni Pei, Fan Hongrui, Pirajno F, Lai Yong, Su Wenchao, Zhang Hui. 2007. Diagnostic fluid inclusions of different types hydrothermal gold deposits. Acta Petrologica Sinica, 23(9): 2085~2108 (in Chinese with English abstract).
-
Chen Zhenyan, Yan Huo, Li Junsheng, Zhang Ge, Zhang Zhanwen, Liu Baozhu. 1999. Relationship between tertiary volcanic rocks and hydrocarbons in the Liaohe basin, People's Republic of China. AAPG Bulletin, 83(6): 1004~1014.
-
Cooper J R, Crelling J C, Rimmer S M, Whittington A G. 2007. Coal metamorphism by igneous intrusion in the Raton basin, CO and NM: implications for generation of volatiles. International Journal of Coal Geology, 71(1): 15~27.
-
Cui Zhuang, Hou Bing, Fu Shihao, Lv Jiaxin. 2022. Fractures cross-layer propagation characteristics of integrated fracturing in shale oil tight reservoir. Fault Block Oil & Gas Field, 29(1): 111~117 (in Chinese with English abstract).
-
Das R, Zhang Y, Schaubs P, Cleary P W. 2014. Modelling rock fracturing caused by magma intrusion using the smoothed particle hydrodynamics method. Computational Geosciences, 18(6): 927~947.
-
Delaney P T. 1982. Rapid intrusion of magma into wet rock: groundwater flow due to pore pressure increases. Journal of Geophysical Research, 87(B9): 7739~7756.
-
Delaney P T. 1987. Heat transfer during emplacement and cooling of mafic dykes. In: Mafik Dike Swarms (Halls H L, Fahring W F, eds). Geologic Association, Canada. Special Paper, 34: 31~46.
-
Duuring P, Hagemann S G, Groves D I. 2000. Structural setting, hydrothermal alteration, and gold mineralisation at the Archaean syenite-hosted Jupiter deposit, Yilgarn Craton, western Australia. Mineralium Deposita, 35(5): 402~421.
-
Engvik A K, Bertram A, Kalthoff J F, Stöckhert B, Austrheim H, Elvevold S. 2005. Magma-driven hydraulic fracturing and infiltration of fluids into the damaged host rock, an example from Dronning Maud Land, Antarctica. Journal of Structural Geology, 27(5): 839~854.
-
Feng Jianwei, Shi Shi, Zhou Zhaohua, Li Xizhe, Luo Peng. 2019. Characterizing the influence of interlayers on the development and distribution of fractures in deep tight sandstones using finite element method. Journal of Structural Geology, 123: 81~95.
-
Feng Yuntian, Owen D R J. 2014. Discrete element modelling of large scale particle systems—I: exact scaling laws. Computational Particle Mechanics, 1(2): 159~168.
-
Fu F Q, McInnes B I A, Evans N J, Davies P J. 2010. Numerical modeling of magmatic-hydrothermal systems constrained by U-Th-Pb-He time-temperature histories. Journal of Geochemical Exploration, 106(1-2): 90~109.
-
Fu Xiaofang, Wang Chenghong, Yuan Linping, Liang Bin, Hao Xuefeng, Pan Meng. 2014. Achievements in the investigation and evaluation of spodumene sesources at Jiajika in Sichuan, China. Geological Survey of China, 1(3): 37~43 (in Chinese with English abstract).
-
Galushkin Y I. 1997. Thermal effects of igneous intrusions on maturity of organic matter: a possible mechanism of intrusion. Organic Geochemistry, 26(11-12): 645~658.
-
Giudicepietro F, Macedonio G, Auria L D, Martini M. 2016. Insight into vent opening probability in volcanic calderas in the light of a sill intrusion model. Pure and Applied Geophysics, 173(5): 1703~1720.
-
Gray G G, Morgan J K, Sanz P F. 2014. Overview of continuum and particle dynamics methods for mechanical modeling of contractional geologic structures. Journal of Structural Geology, 59: 19~36.
-
Hayba D O, Ingebritsen S E. 1997. Multiphase groundwater flow near cooling plutons. Journal of Geophysical Research: Solid Earth, 102(B6): 12235~12252.
-
Hobbs B E, Ord A. 2010. The mechanics of granitoid systems and maximum entropy production rates. Philosophical Transactions of the Royal Society a Mathematical Physical and Engineering Sciences, 368(1910): 53~93.
-
Hou Zengqian, Mo Xuanxue, Gao Yongfeng, Qu Xiaoming, Meng Xiangjin. 2003. Adakite, a possible host rock for porphyry copper deposits: case studies of porphyry copper belts in Tibetan Plateau and in Northern Chile. Mineral Deposits, 22(1): 1~12 (in Chinese with English abstract).
-
Hu Ruizhong, Zhou Meifu. 2012. Multiple Mesozoic mineralization events in South China—an introduction to the thematic issue. Mineralium Deposita, 47(6): 579~588.
-
Jaeger J C. 1959. Temperatures outside a cooling intrusive sheet. American Journal of Science, 257(1): 44~54.
-
Krynauw J R, Behr H J, Van Den Kerkhof A M. 1994. Sill emplacement in wet sediments: fluid inclusion and cathodoluminescence studies at Grunehogna, western Dronning Maud Land, Antarctica. Journal of the Geological Society, 151(5): 777~794.
-
Li Changsheng, Yin Hongwei, Wu Chuang, Zhang Yingchun, Zhang Jiaxing, Wu Zhenyun, Wang Wei, Jia Dong, Guan Shuwei, Ren Rong. 2021. Calibration of the discrete element method and modeling of shortening experiments. Frontiers in Earth Science, 9: 1~12.
-
Li Xia, Wang Qin, Zhang Wen, Yin Hongwei. 2016. Contact metamorphism of shales intruded by a granite dike: implications for shale gas preservation. International Journal of Coal Geology, 159: 96~106.
-
Liu Chao, Xie Qingbin, Wang Guiwen, Cui Yu, Zhang Chujun. 2015. The influence of igneous intrusion to detrital reservoir: advances and outlook. Advance in Earth Sciences, 30(6): 654~667 (in Chinese with English abstract).
-
Liu Chun, Xu Qiang, Shi Bin, Deng Shang, Zhu Honghu. 2017. Mechanical properties and energy conversion of 3D close-packed lattice model for brittle rocks. Computers & Geosciences, 103: 12~20.
-
Liu Chun, Le Tiancheng, Shi Bin, Zhu Yao. 2020. Discussion on three major problems of engineering application of the particle discrete element method. Chinese Journal of Rock Mechanics and Engineering, 39(6): 1142~1152 (in Chinese with English abstract).
-
Liu Honglin, Wang Hongyan, Liu Renhe, Zhao Qun, Lin Yingji. 2010. China shale gas resources and prospect potential. Acta Geologica Sinica, 84(9): 1374~1378 (in Chinese with English abstract).
-
Lovering T S. 1935. Theory of heat conduction applied to geological problems. Geological Society of America Bulletin, 46(1): 69~94.
-
Ma Yemu, Lu Xiancai, Zhang Xuefen, Li Xiaozhao, Hu Wenxuan, Wang Kai. 2013. A numerical simulation of the heat transfer in granite intrusion-mudstone contact zone and its geological implication: a case study from eastern Guangdong Province, China. Geological Journal of China Universities, 19(2): 307~315 (in Chinese with English abstract).
-
Malfait W J, Seifert R, Petitgirard S, Perrillat J P, Mezouar M, Ota T, Nakamura E, Lerch P, Sanchez-Valle C. 2014. Supervolcano eruptions driven by melt buoyancy in large silicic magma chambers. Nature Geoscience, 7(2): 122~125.
-
Mao Jingwen, Cheng Yanbo, Chen Maohong, Franco P. 2012. Major types and time-space distribution of Mesozoic ore deposits in South China and their geodynamic settings. Mineralium Deposita, 48(3): 267~294.
-
Ni Pei, Chi Zhe, Pan Junyi. 2020. An integrated investigation of ore-forming fluid evolution in porphyry and epithermal deposits and their implication on exploration. Earth Science Frontiers, 27(2): 60~78 (in Chinese with English abstract).
-
Nie Haikuan, He Zhiliang, Liu Guangxiang, Zhang Guangrong, Lu Zhiyuan, Li Donghui, Sun Chuanxiang. 2020. Status and direction of shale gas exploration and development in China. Journal of China University of Mining & Technology, 49(1): 13~35 (in Chinese with English abstract).
-
Othman R, Arouri K R, Ward C R, McKirdy D M. 2001. Oil generation by igneous intrusions in the northern Gunnedah basin, Australia. Organic Geochemistry, 32(10): 1219~1232.
-
Ouyang Yongpeng, Rao Jianfeng, Yao Zaiyu, Zhou Xianrong, Chen Guohua. 2018. Mineralization and prospecting direction of the Zhuxi type skarn deposit. Geological Science and Techology Information, 37(3): 148~158 (in Chinese with English abstract).
-
Pan Jiping, Qiao Dewu, Li Shizheng, Zhou Dongsheng, Xu Linfeng, Zhang Mengyin, Song Xiuyan. 2011. Shale-gas geological conditions and exploration prospect of the Paleozoic marine strata in lower Yangtze area, China. Geological Bulletin of China, 30(2-3): 337~343 (in Chinese with English abstract).
-
Peacock S M. 1990. Numerical simulation of regional and contact metamorphism using the Macintosh microcomputer. Journal of Geology, 38: 132~137.
-
Qin Kezhang, Zhou Qifeng, Zhao Junxing, He Changtong, Liu Xiaochi, Shi Ruizhe, Liu Yuchao. 2021a. Be-rich mineralization features of Himalayan leucogranite belt and prospects for lithium-bearing pegmatites in higher altitudes. Acta Geologica Sinica, 95(10): 3039~3053 (in Chinese with English abstract).
-
Qin Kezhang, Zhou Qifeng, Tang Dongmei, Wang Chunlong, Zhu Liqun. 2021b. The emplacement mechanism, melt-fluid evolution, rare-element metallogenesis and puzzles of the Koktokay No. 3 pegmatite rare elemental deposit, Altai. Acta Geologica Sinica, 95(10): 3146~3162 (in Chinese with English abstract).
-
Qiu Longwei, Xi Qingfu, Liu Kuiyuan. 2000. A study on attitude of insidious diabase and its lithoface zone in Luo151 well field in Zhanhua depression, Shandong. Acta Petrologica Sinica, 16(3): 413~419 (in Chinese with English abstract).
-
Raymond A C, Murchison D G. 1988. Development of organic maturation in the thermal aureoles of sills and its relation to sediment compaction. Fuel, 67(12): 1599~1608.
-
Richards J P. 2015. Tectonic, magmatic, and metallogenic evolution of the Tethyan orogen: from subduction to collision. Ore Geology Reviews, 70: 323~345.
-
Rubin A M. 1995. Propagation of magma-filled cracks. Annual Review of Earth and Planetary Sciences, 23(1): 287~336.
-
Santos R V, Dantas E L, deOliveira C G, deAlvarenga C J S, Anjos C W D d, Guimarães E M, Oliveira F B. 2009. Geochemical and thermal effects of a basic sill on black shales and limestones of the Permian Irati Formation. Journal of South American Earth Sciences, 28(1): 14~24.
-
Simoneit B R T, Brenner S, Peters K E, Kaplan I R. 1981. Thermal alteration of Cretaceous black shale by diabase intrusions in the Eastern Atlantic—II. effects on bitumen and kerogen. Geochimica et Cosmochimica Acta, 45(9): 1581~1602.
-
Svensen H, Jamtveit B. 1998. Contact metamorphism of shales and limestones from the Grua area, the Oslo Rift, Norway: a phase-petrological study. Norsk Geologisk Tidsskrift, 78: 81~98.
-
Svensen H, Planke S, Chevallier L, Malthe-Sørenssen A, Corfu F, Jamtveit B. 2007. Hydrothermal venting of greenhouse gases triggering Early Jurassic global warming. Earth and Planetary Science Letters, 256(3-4): 554~566.
-
Wang Dayong, Lu Xiancai, Zhang Xuejun, Xu Shijin, Hu Wenxuan, Wang Liangshu. 2007. Heat model analysis of wall rocks below a diabase sill in Huimin sag, China compared with thermal alteration of mudstone to carbargilite and hornfels and with increase of vitrinite reflectance. Geophysical Research Letters, 34(16): https: //doi. org/10. 1029/2007GL030314.
-
Wang Dayong, Lu Xiancai, Song Yongchen, Shao Rong, Qi Tian. 2010. Influence of the temperature dependence of thermal parameters of heat conduction models on the reconstruction of thermal history of igneous-intrusion-bearing basins. Computers & Geosciences, 36(10): 1339~1344.
-
Wang Dayong, Lu Xiancai, Xu Shijin, Hu Wenxuan, Qi Tian. 2011. Heat-transfer-model analysis of thermal effect of intrusive sills on organic-rich host rocks in sedimentary basins. Journal of Nanjing University, 47(1): 45~50 (in Chinese with English abstract).
-
Wang Dayong, Song Yongchen, Liu Weiguo, Zhao Minglong, Qi Tian. 2011. Numerical investigation of the effect of volatilization and the supercritical state of pore water on maturation of organic matter in the vicinity of igneous intrusions. International Journal of Coal Geology, 87(1): 33~40.
-
Wang Dayong, Manga M. 2015. Organic matter maturation in the contact aureole of an igneous sill as a tracer of hydrothermal convection. Journal of Geophysical Research: Solid Earth, 120(6): 4102~4112.
-
Wang Guoguang, Ni Pei, Pan Junyi. 2020. Fluid characteristics of granite-related ore forming systems. Bulletin of Mineralogy Petrology and Geochemistry, 39(3): 463~471 (in Chinese with English abstract).
-
Wang He, Li Pei, Ma Huadong, Zhu Binyu, Qiu Lin, Zhang Xiaoyu, Dong Rui, Zhou Kailin, Wang Min, Wang Qian, Yan Qinghe, Wei Xiaopeng, He Bin, Lu He, Gao Hao. 2017. Discovery of the Bailongshan superlarge lithium-rubidium deposit in Karakorum, Hetian, Xinjiang, and its prospecting implication. Geotectonica et Metallogenia, 41(6): 1053~1062 (in Chinese with English abstract).
-
Wang Man, Wang Yingwei, Xue Linfu. 2012. Finite element simulation of effect of magmatic intrusion on maturity of organic matter in surrounding rock. Fault-Block Oil & Gas Field, 19(2): 172~176 (in Chinese with English abstract).
-
Wang Min, Lu Shuangfang, Xue Haitao, Wu Jing, Liu Dawei. 2010. The effects of magmatic intrusions on the maturation of organic matter and its numerical simulation. Acta Petrologica Sinica, 26(1): 177~184 (in Chinese with English abstract).
-
Wang Rucheng, Wu Fuyuan, Xie Lei, Liu Xiaochi, Wang Jiamin, Yang Lei, Lai Wen, Liu Chen. 2017. A preliminary study of rare-metal mineralization in the Himalayan leucigranite belts, South Tibet. Science China: Earth Sciences, 47(8): 871~880 (in Chinese with English abstract).
-
Wang Tao, Wang Xiaxia, Li Wuping. 2000. Evaluation of multiple emplacement mechanisms: the Huichizi granite pluton, Qinling orogenic belt, central China. Journal of Structural Geology, 22(4): 505~518.
-
Wang Yue, Liu Chun, Liu Xiaolei, Liu Hui, Li Yasha. 2021. Discrete element numerical simulation of the accumulation process of wave-induced pore water pressure in the seabed. Acta Oceanologica Sinica, 43(11): 88~95 (in Chinese with English abstract).
-
Winter J D. 2001. An Introduction to Igneous and Metamorphic Petrology. New Jersey: Prentice-Hall Inc, 1~697.
-
Wu Fuyuan, Wang Rucheng, Liu Xiaochi, Xie Lei. 2021. New breakthroughs in the studies of Himalayan rare-metal mineralization. Acta Petrologica Sinica, 37(11): 3261~3276 (in Chinese with English abstract).
-
Xu Zhiqin, Wang Rucheng, Zhao Zhongbao, Fu Xiaofang. 2018. On the structural backgrounds of the large-scale hard-rock type lithium ore belts in China. Acta Geologica Sinica, 92(6): 1091~1106 (in Chinese with English abstract).
-
Yang Zongfeng, Luo Zhaohua, Lu Xinxiang, Huang Fan, Chen Bihe. 2011. Geological characteristics and prospecting potential of the Shanggusi porphyry molybdenum deposit in the east Qinling. Geology and Prospecting, 47(6): 1077~1090 (in Chinese with English abstract).
-
Zhang Hongyong, Zhang Yanjie, Liu Chun, Shi Bin, Cao Zheng. 2021. Analysis of grouting treatment for convergence deformation of meto tunnel based on discrete element pore density flow method. Hazard Control in Tunnelling and Underground Engineering, 3(3): 100~110 (in Chinese with English abstract).
-
Zhang Jian, Shi Yaolin. 1997. Thermal simulation of magmatic intrusion in sedimentary basin. Progress in Geophysics, 12(3): 55~64 (in Chinese with English abstract).
-
Zhang Wen, Wang Qin, Ye Jianguo, Zhou Jie. 2017. Fracture development and fluid pathways in shales during granite intrusion. International Journal of Coal Geology, 183: 25~37.
-
Zhang Wen, Wang Qin. 2018. Permeability anisotropy and gas slippage of shales from the Sichuan basin in South China. International Journal of Coal Geology, 194: 22~32.
-
Zhao Chongbin, Hobbs B E, Ord A, Peng Shenglin. 2008. Particle simulation of spontaneous crack generation associated with the laccolithic type of magma intrusion processes. International Journal for Numerical Methods in Engineering, 75(10): 1172~1193.
-
Zhao Junxing, Qin Kezhang, Li Guangming, Li Jinxiang, Xiao Bo, Chen Lei. 2012. Geochemistry and petrogenesis of granitoids at Sharangeocene porphyry Mo deposit in the main-stage of India-Asia continental collision, Northern Gangdese, Tibet. Resource Geology, 62(1): 84~98.
-
Zhou Zongqing, Ranjith P G, Yang Weimin, Shi Shaoshuai, Wei C C, Li Z H. 2019. A new set of scaling relationships for DEM-CFD simulations of fluid-solid coupling problems in saturated and cohesiveless granular soils. Computational Particle Mechanics, 6(4): 657~669.
-
蔡周荣, 夏斌, 黄强太, 万志峰. 2015. 上、下扬子区古生界页岩气形成和保存的构造背景对比分析. 天然气地球科学, 26(8): 1446~1454.
-
陈骏, 王汝成, 朱金初, 陆建军, 马东升. 2013. 南岭多时代花岗岩的钨锡成矿作用. 中国科学: 地球科学, 44(1): 111~121.
-
陈衍景. 2010. 秦岭印支期构造背景、岩浆活动及成矿作用. 中国地质, 37(4): 854~865.
-
陈衍景, 倪培, 范宏瑞, F Pirajno, 赖勇, 苏文超, 张辉. 2007. 不同类型热液金矿系统的流体包裹体特征. 岩石学报, 23(9): 2085~2108.
-
崔壮, 侯冰, 付世豪, 吕嘉昕. 2022. 页岩油致密储层一体化压裂裂缝穿层扩展特征. 断块油气田, 29(1): 111~117.
-
付小方, 侯立玮, 王登红, 袁蔺平, 梁斌, 郝雪峰, 潘蒙. 2014. 四川甘孜甲基卡锂辉石矿矿产调查评价成果. 中国地质调查, 1(3): 37~43.
-
侯增谦, 莫宣学, 高永丰, 曲晓明, 孟祥金. 2003. 埃达克岩: 斑岩铜矿的一种可能的重要含矿母岩——以西藏和智利斑岩铜矿为例. 矿床地质, 22(1): 1~12.
-
刘超, 谢庆宾, 王贵文, 崔宇, 张楚珺. 2015. 岩浆侵入作用影响碎屑围岩储层的研究进展与展望. 地球科学进展, 30(6): 654~667.
-
刘春, 乐天呈, 施斌, 朱遥. 2020. 颗粒离散元法工程应用的三大问题探讨. 岩石力学与工程学报, 39(6): 1142~1152.
-
刘洪林, 王红岩, 刘人和, 赵群, 林英姬. 2010. 中国页岩气资源及其勘探潜力分析. 地质学报, 84(9): 1374~1378.
-
马野牧, 陆现彩, 张雪芬, 李晓昭, 胡文瑄, 汪恺. 2013. 花岗岩侵入体-泥质围岩热传输过程的数值模拟及其地质意义: 以粤东典型接触带剖面为例. 高校地质学报, 19(2): 307~315.
-
倪培, 迟哲, 潘君屹. 2020. 斑岩型和浅成低温热液型矿床成矿流体与找矿预测研究: 以华南若干典型矿床为例. 地学前缘, 27(2): 60~78.
-
聂海宽, 何治亮, 刘光祥, 张光荣, 卢志远, 李东晖, 孙川翔. 2020. 中国页岩气勘探开发现状与优选方向. 中国矿业大学学报, 49(1): 13~35.
-
欧阳永棚, 饶建锋, 尧在雨, 周显荣, 陈国华. 2018. 朱溪式矽卡岩型矿床成矿作用及找矿方向. 地质科技情报, 37(3): 148~158.
-
潘继平, 乔德武, 李世臻, 周东升, 许林峰, 张梦颖, 宋修艳. 2011. 下扬子地区古生界页岩气地质条件与勘探前景. 地质通报, 30(Z1): 337~343.
-
秦克章, 周起凤, 唐冬梅, 王春龙, 朱丽群. 2021a. 阿尔泰可可托海3号脉花岗伟晶岩侵位机制、熔-流体演化、稀有金属富集机理及待解之谜. 地质学报, 95(10): 3039~3053.
-
秦克章, 周起凤, 赵俊兴, 何畅通, 刘小驰, 施睿哲, 刘宇超. 2021b. 喜马拉雅淡色花岗岩带伟晶岩的富铍成矿特点及向更高处找锂. 地质学报, 95(10): 3146~3162.
-
邱隆伟, 席庆福, 刘魁元. 2000. 山东沾化凹陷罗151井区隐伏辉绿岩的产状及岩相带划分. 岩石学报, 16(3): 413~419.
-
王大勇, 陆现彩, 徐士进, 胡文瑄, 齐天. 2011a. 沉积盆地内侵入岩席对富含有机质围岩热影响的热传输模型研究. 南京大学学报(自然科学版), 47(1): 45~50.
-
王国光, 倪培, 潘君屹. 2020. 花岗质岩石相关成矿系统的流体作用. 矿物岩石地球化学通报, 39(3): 463~471.
-
王核, 李沛, 马华东, 朱炳玉, 邱林, 张晓宇, 董瑞, 周楷麟, 王敏, 王茜, 闫庆贺, 魏小鹏, 何斌, 卢鸿, 高昊. 2017. 新疆和田县白龙山超大型伟晶岩型锂铷多金属矿床的发现及其意义. 大地构造与成矿学, 41(6): 1053~1062.
-
王满, 王英伟, 薛林福. 2012. 岩浆侵入对围岩中有机质成熟度影响的有限元模拟. 断块油气田, 19(2): 172~176.
-
王民, 卢双舫, 薛海涛, 武静, 刘大为. 2010. 岩浆侵入体对有机质生烃(成熟)作用的影响及数值模拟. 岩石学报, 26(1): 177~184.
-
王汝成, 吴福元, 谢磊, 刘小驰, 王佳敏, 杨雷, 赖文, 刘晨. 2017. 藏南喜马拉雅淡色花岗岩稀有金属成矿作用初步研究. 中国科学: 地球科学, 47(8): 871~880.
-
王岳, 刘春, 刘晓磊, 刘辉, 李亚沙. 2021. 波浪作用下海床孔压累积过程离散元数值模拟. 海洋学报, 43(11): 88~95.
-
吴福元, 王汝成, 刘小驰, 谢磊. 2021. 喜马拉雅稀有金属成矿作用研究的新突破. 岩石学报, 37(11): 3261~3276.
-
许志琴, 王汝成, 赵中宝, 付小方. 2018. 试论中国大陆“硬岩型”大型锂矿带的构造背景. 地质学报, 92(6): 1091~1106.
-
杨宗锋, 罗照华, 卢欣祥, 黄凡, 陈必河. 2011. 东秦岭尚古寺斑岩钼矿地质特征及成矿潜力分析. 地质与勘探, 47(6): 1077~1090.
-
张鸿勇, 张艳杰, 刘春, 施斌, 曹政. 2021. 基于离散元孔隙密度流法的地铁隧道收敛变形注浆整治分析. 隧道与地下工程灾害防治, 3(3): 100~110.
-
张健, 石耀霖. 1997. 沉积盆地岩浆侵入的热模拟. 地球物理学进展, 12(3): 55~64.
-
摘要
岩浆侵位会引起围岩的接触变质作用和变形,虽然围岩中的破裂提供了流体活动和元素迁移的重要通道,但是大多数岩浆侵位模型忽略了围岩中热对流的影响。本文总结了岩浆侵位的数值模拟原理和常用模拟方法,然后使用离散元软件MatDEM建立了二维的双层围岩模型和均质围岩模型,用孔隙密度流法模拟花岗质岩浆侵位和冷却过程中岩浆与围岩之间的流-固-热-力耦合过程。结果表明岩浆侵位过程中围岩的裂隙发育和热对流对接触变质晕具有重要影响,围岩的裂隙发育和热传输模式可分为三个阶段:① 岩浆侵位初期的挤压力使围岩中产生广泛分布的、呈径向展布的剪裂隙,围岩以热传导和孔隙渗流为主导;② 在持续的孔隙流体压下,径向裂隙连通形成主干张性断裂并向上扩展,成为熔/流体迁移的重要通道以及伟晶岩型和热液型矿床的成矿空间,通道流和局部热对流控制了围岩的热传输;③ 在岩浆侵位后期,岩浆房附近围岩中的孔隙流体压增大,在侵入体与围岩的接触带形成大量张裂,加强了变质晕内的热对流,有助于矽卡岩型矿床的形成。与只有热传导的模型对比,热对流使围岩中变质晕的宽度减小。接触变质晕的几何形态受侵入体的形态控制,但是变质晕宽度在空间上有显著差异。本研究为重建岩浆侵位过程中的变质-变形-成矿作用提供了新方法。
Abstract
Emplacement of granitic magma causes contact metamorphism and brittle deformation in host rocks. These fractures provide important pathways for fluid activity and element migration. However, most magma emplacement models ignored the influence of thermal convection. Here we summarize numerical simulation principles of magma emplacement and compare different simulation methods. Using the MatDEM discrete element software, we establish a double-layer host rock model and a homogeneous host rock model. The pore-density flow method of MatDEM was applied to simulate the fluid-solid-heat-stress coupling process during intrusion and cooling of granitic magma. The results reveal three stages of the fracture development and heat transfer patterns during magma emplacement. ① During initial intrusion, compression of magma produces the widespread, radial distribution of shear fractures in host rocks. Heat conduction and infiltration are dominant in host rocks. ② Under continuous pore fluid pressure, the radial shear fractures are connected and form major extensional fractures, which become important pathways of melt and fluids and provide space for pegmatitic and hydrothermal ore deposits. Channel flow and localized thermal convection control heat transportation in host rocks. ③ In the last stage, increase of the pore fluid pressure in host rocks around the magma chamber triggers large amounts of tensile fractures, which enhance the thermal convection in the contact aureole and promoting skarn-type mineralization. Compared with the thermal conduction model, initial thermal convection will decrease the width of the contact aureole. The geometric shape of a contact aureole is controlled by the shape of intrusion, but its width shows remarkable spatial variations. Our modeling results provide a new way to reconstruct metamorphism, deformation and mineralization during magma emplacement.
