石英脉型黑钨矿床是重要的钨矿床类型之一，主要分布于南岭成矿带的赣南、湘南、粤北、桂北等地^[1-2]。钨矿脉主要分布在碱长花岗岩顶部接触带附近^[3]。大多数钨矿脉近垂直产出，垂向延伸可达1000m^[2]。石英型黑钨矿床矿脉围岩蚀变呈带状分布在矿脉两侧，围岩蚀变宽度随深度递减，且这种变化趋势不受岩性影响^{[1, 4]}。多位学者认为，形成黑钨矿石英脉的成矿流体不是纯热水溶液，而是富含SiO₂和挥发分的岩浆热液过渡性流体^[4-7]。这种岩浆热液过渡性流体比水溶液具有更高的粘度和密度。形成钨矿脉的过渡性流体在深部偏向熔体，在浅部偏向热液，因而蚀变宽度上宽下窄。Audétat等^[8]的实验表明，只有当熔体的摩尔分数很高时，含水熔体的粘度才会显著高于纯水的粘度。然而，除常海亮等^[9]在西华山钨矿脉的绿柱石中发现熔体包裹体外，大多数矿石矿物和脉石矿物的流体包裹体研究结果表明，钨成矿流体是NaCl-H₂O±CO₂体系^[10-22]。因而，形成上述蚀变特征可能另有他因。岩浆热液矿床的形成涉及复杂的化学反应和热液运移过程^[23-25]。前人研究多偏重于前者，对热液运移过程的研究较少^[26-28]。本文利用有限元方法模拟热液运移和化学组分从裂隙向邻近围岩扩散过程，发现围岩孔隙度随深度增加而递减可能是钨矿床形成上述蚀变特征的有效机制。

1.   石英脉型钨矿床地质和地球化学特征

中国石英脉型钨矿床主要分布于南岭成矿带的赣南、湘南、粤北、桂北等地（图 1）。南岭地区大规模钨成矿作用与燕山早期的花岗岩岩浆活动有关，成矿年龄大多集中在中晚侏罗世(165~ 150Ma) ^[29]。该类矿床常产在碱长花岗岩顶上带附近^[3]。主要的矿石矿物为黑钨矿、白钨矿、锡石等，脉石矿物有石英、长石、萤石、黄铁矿等^[1]。

图  1  南岭成矿带钨矿分布示意图^[29]

Figure  1.  Distribution of tungsten deposits in the Nanling Mountains, South China

下载: 全尺寸图片 幻灯片

钨成矿流体主要属于中高温、中低盐度的NaCl-H₂O±CO₂体系^[30]。首先，成矿流体成分以H₂O为主，CO₂含量低，少量CO、CH₄、N₂、H₂。在主成矿阶段成矿流体主要来自岩浆水，在成矿晚期成矿流体中混有大量大气降水。其次，寄主矿物（石英、黑钨矿、绿柱石、锡石等）捕获的流体包裹体在常温下以富液相的气液两相包裹体为主，纯气相和含CO₂的三相包裹体少见，表明溶液不曾沸腾或局部沸腾，压力可达90~160MPa。流体包裹体的均一温度一般为160~390℃，矿石矿物（黑钨矿、锡石等）的均一温度比石英的高40~100℃。最后，成矿流体盐度较低, 成矿早期流体盐度一般为7%~10%NaCl，晚期逐渐降低至1%^{[1-2, 10-22, 31-35]}。

钨矿床围岩类型以寒武系和奥陶系居多，岩性以沉积(变质)碎屑岩为主^{[2, 36]}。蚀变类型包括云英岩化、硅化、黑云母化、电气石化等^[2]。含钨石英脉的围岩蚀变有以下特点^[1]：①围岩蚀变呈带状分布在矿脉两侧（图 2），蚀变带宽几毫米至几十厘米，带状分布特点不受岩性影响；②围岩蚀变宽度上宽下窄，蚀变宽度变化趋势不受岩性影响；③矿脉以石英为主，但石英脉体两侧的硅化较弱，且上部硅化比下部强。

图  2  南岭地区石英脉型钨矿典型矿脉

a—江西大吉山钨矿深部脉体，标高417m，脉体走向为近东西向，陡倾，向上拍摄；b—江西漂塘钨矿深部大脉，标高268m，脉体走向为近东西向，陡倾，向上拍摄。Qtz—石英；Wol—黑钨矿

Figure  2.  Typical veins at deeper levels of tungsten deposits in the Nanling Mountains

下载: 全尺寸图片 幻灯片

反应性流体（如成矿热液）流经岩石裂隙或孔隙时，会改变邻近围岩的化学成分，形成蚀变晕^[37]。蚀变晕的宽度与岩石、流体的物理化学性质、反应速率、反应时间等参数密切相关^[37-39]。本文利用热液运移数值模拟研究化学组分从裂隙向邻近围岩扩散过程，为解释石英脉型钨矿床蚀变特征提供新思路。

2.   热液运移数值模拟

2.1   数学理论

多孔介质流体流动和热传输用流体连续性方程、达西定律和热传输方程定量描述，这些偏微分方程组用邢会林教授开发的PANDAS (Parallel Adaptive Nonlinear Deformation Analysis System)来求解。PANDAS是一套基于有限元方法(FEM)和格子Boltzmann方法(LBM)开发的创新性软件平台，已经成功用于模拟地下含裂隙非均质孔隙材料/地层中“应力变形/破坏-流体流动-热传导-化学反应”等多物理场高度非线性耦合的问题，解决从微孔隙尺度到实验室及油气田尺度，甚至全球等多尺度的地球科学和地质资源工程中的主要科学问题和工程技术挑战^[40-49]。

NaCl水溶液的密度和粘度与温度、压力和盐度有关。Batzle等提出的盐水密度模型在热液条件下精度较高^[50]，因而本文的模拟计算使用了该密度模型。笔者提出一个适于大规模数值计算的盐水粘度模型^[43]，用于计算热液的粘度。

某一组分在多孔介质的有效扩散速率D_e与该组分在水溶液的扩散速率D_s和孔隙度ϕ有关^[51]：

扩散速率D_s服从Stokes-Einstein方程^[52-53]。

Cathles等^[37]提出垂直于脉体的蚀变晕宽度计算公式：

其中，Z为时间t内形成的蚀变晕宽度，G是一个无量纲数，其值为改变1m³岩石所需要的流体体积。本文利用Péclet数（Pe）表征流体传输机制^[54]，当Pe < 1时，流体流动的扩散比例较高，流体内的化学反应以反应-扩散过程为主；当Pe > 1时，流体流动的对流比例升高，化学反应以反应-扩散-对流过程为主^[55]。

2.2   热液运移模型

本文主要根据以下背景资料建立热液运移模型：①流体包裹体测温研究表明，钨成矿流体的压力较高，可达到90~160MPa，平均值为125MPa；②钨成矿热液来自深部含水岩浆房，岩浆顶上带与围岩接触部位发育大量的裂隙，这些裂隙为热液运移提供构造通道和成矿空间^[56-58]。

根据石英脉型钨矿床的地质地球化学特征，本文建立了二维热液运移模型（图 3）。Z轴代表重力方向，X轴为垂直脉带的水平方向。为了减少网格数量和计算时间，本文模型只截取深度3.8~5.0km的部分，模型尺寸为1.2km×1.2km。模型包含裂隙区和围岩2个单元，石英脉型钨矿床脉带宽度从几十米至几百米不等^[2]，模型中裂隙区宽度为20m。选取2条剖面用于研究裂隙带及其紧邻围岩的热液运移，其中一条剖面在4.9km，另一条剖面在4.1km处。

图  3  石英脉型钨矿床热液运移数值模型

Figure  3.  Numerical model of hydrothermal flow at the tungsten deposits in the Nangling Mountains

下载: 全尺寸图片 幻灯片

设定模型底部边界条件为压力125MPa和温度400℃，顶部边界条件为压力38MPa和温度114℃。根据现今地热储库观测数据^[59]，设定初始温度梯度为30℃/km，初始压力梯度为10MPa/km，模型中热液盐度为10% NaCl。硅化是南岭石英脉钨矿床较普遍的蚀变类型之一，本文将硅从裂隙到邻近围岩的扩散作为研究对象。硅在25℃稀溶液的扩散系数约为1.0× 10^-9m²/s^[60]。

本文设计2个数值实验，2个实验中裂隙带比围岩拥有更高的孔隙度和渗透率。在实验1中，围岩和裂隙带的孔隙度和渗透率均一分布，裂隙带孔隙度为10%，渗透率为1.0 × 10^-12m²；围岩孔隙度为0.01，渗透率为1.0× 10^-16m²。在实验2中，围岩的孔隙度和渗透率与深度呈负指数关系^[61]，渗透率在模型截取的深度范围内变化约一个数量级。

3.   模拟结果

3.1   数值实验1

在数值实验1中，围岩的孔隙度和渗透率均一分布。4.9km剖面和4.1km剖面的初始温度分别为147 ℃和123℃。热液运移1a后，裂隙带的温度显著高于邻近围岩的温度，但远离裂隙带的围岩温度变化较小（图 4），在4.9km剖面上，紧邻裂隙带的围岩温度达到340℃；在4.1km剖面上，紧邻裂隙带的围岩温度为320℃。在2条剖面中，裂隙带内的Péclet数为1.7~1.8，围岩的Péclet数远低于1。

图  4  数值实验1温度(a)和Péclet数(b)的分布

Figure  4.  The temperature (a) and Péclet number (b) in the first numerical experiment

下载: 全尺寸图片 幻灯片

由于4.9km剖面的初始温度高于4.1km剖面的初始温度，前者的初始有效扩散系数高于后者。在紧邻裂隙的围岩内，硅有效扩散系数在热液运移1a后有显著提升：在4.9km剖面最大值达2.3×10^-13m²/s，在4.1km剖面最大值为2.1×10^-13m²/s（图 5）。热液运移2a后，围岩在4.1km剖面的有效扩散系数仍然低于围岩在4.9km的有效扩散系数。

图  5  数值实验1热液运移1a后(a)和2a后(b)硅有效扩散系数在裂隙带附近的分布

Figure  5.  The effective diffusion coefficients of silica in the wallrock adjacent to fracture zones after one year (a) and two years (b) in the first numerical experiment

下载: 全尺寸图片 幻灯片

3.2   数值实验2

在数值实验2中，围岩的孔隙度和渗透率随深度增加而按指数规律递减，其他参数与数值实验1一致。在该实验中，围岩孔隙度从3.8km的1.17%降至5.0km的0.58%，对应的渗透率由10^-15.8m²降至10^-16.7m²。围岩在4.9km的孔隙度为0.64%，渗透率为10^-16.6m²。围岩在4.1km的孔隙度为0.99%，渗透率为10^-16.0m²。

热液运移1a后，紧邻裂隙带的围岩在4.9km的温度达到315℃，而相应的温度在4.1km为230℃（图 6）。显著低于数值实验1中紧邻裂隙带的围岩温度，裂隙带的Péclet数为1.8~2.0，围岩的Péclet数低于1，说明实验2中热液运移的流体传输机制与实验1相似。

图  6  数值实验2温度(a)和Péclet数(b)的分布

Figure  6.  The temperature (a) and Péclet number (b) in the second numerical experiment

下载: 全尺寸图片 幻灯片

该实验的初始有效扩散系数与实验1有明显区别。围岩在4.9km的初始有效扩散系数为3.2 × 10^-14m²/s，约为4.1km剖面的初始有效扩散系数(6.0 × 10^-14m²/s)的一半。热液运移1a后，围岩在4.1km的有效扩散系数最大值为1.4×10^-13m²/s，是4.9km剖面对应值的1.7倍（图 7）。热液运移2a后，围岩在4.1km的有效扩散系数最大值升高至2.2× 10^-13m²/s，是4.9km剖面对应值的2.4倍。

图  7  数值实验2中热液运移1a后(a)和2a后(b)硅有效扩散系数在裂隙带邻近围岩的分布

Figure  7.  The effective diffusion coefficients of silica in the wallrock adjacent to fracture zones after one year (a) and two years (b) in the second numerical experiment

下载: 全尺寸图片 幻灯片

4.   讨论

前人在20世纪80年代发现南岭地区石英脉型钨矿床蚀变宽度呈上宽下窄的特征^[1]。本文利用方程（2）对比不同深度下的蚀变宽度，讨论蚀变宽度上宽下窄的形成机制。由于成矿流体主要来自深部岩浆，假设方程（2）的参数G是一个常数，在给定时间内垂直脉体的硅化宽度与硅的有效扩散系数呈正比例关系。

在数值实验1中，热液运移在紧邻裂隙带的围岩内以扩散为主，反应-扩散过程在化学反应中占主导地位。由于围岩的孔隙度和渗透率是均一分布的常数，硅在围岩的有效扩散系数只决定于硅在水溶液的扩散系数（见公式1）。围岩在深部的温度高于浅部，因而硅在深部围岩中有更高的有效扩散系数，深部易形成较宽的硅化，这显然与石英脉型钨矿床蚀变特征相反。事实上，前人的数值模拟也得到类似的结论：蚀变宽度随着热液运移远离流体源头而递减^{[37, 39]}。然而，这与南岭地区石英脉型钨矿和某些斑岩型铜矿的蚀变特征不符^[38]。通过数值实验1，认为石英脉型钨矿在1km垂直范围内围岩孔隙度可能不是均一分布的。

在数值实验2中，围岩孔隙度和渗透率随深度增大呈指数递减。尽管深部围岩的温度依然高于浅部，由于硅有效扩散系数与孔隙度的平方成正比，随深度递减的孔隙度抵消了温度对于硅有效扩散系数的影响，因而浅部围岩中硅的有效扩散系数高于深部。在这种情况下，裂隙带周围可能形成与南岭石英脉型钨矿床类似的蚀变特征。

本文数值模拟实验的运行时间远小于岩浆热液系统演化的时间尺度，后者可能在数个百万年甚至几十个百万年的时间尺度保持活跃状态^[62-63]。热液矿床成矿动力学过程涉及多时空尺度相互作用的热液流动、岩石变形和化学反应。受岩石变形、化学反应等过程的影响，与矿床形成相关的高渗透率具有局部性、短暂性、周期性等特征^[64-66]。多方面证据表明，上地壳高渗透会在1~1000a内发生显著衰减^[67]。故本文的模拟实验可看作成矿裂隙渗透率近似不变时部分热液充填过程。另外，在本文模型边界条件的控制下，延长模拟实验运行时间会提高围岩整体的温度和硅的扩散速率，但不会改变这2个变量在深部与浅部的相对大小。故数值模拟实验运行时间虽短，但对理解钨矿床蚀变晕形成机制仍具有一定的启发。

石英脉型钨矿床的围岩岩性可能有利于形成随深度递减的孔隙度和渗透率，南岭地区石英脉型钨矿床围岩岩性以沉积(变质)碎屑岩为主^{[2, 29]}。这些围岩较致密，孔隙度也较低。实验表明，由于微小裂隙的闭合，低孔隙度结晶岩石的渗透率随着压力的增加而呈指数规律减小；而高孔隙度沉积岩对压力的敏感性较弱^[68-69]。因而，石英脉型钨矿床围岩岩性可能有利于形成随深度递减的孔隙度和渗透率。野外证据表明，石英脉型钨矿床围岩的裂隙发育程度随深度增加而递减^[56-57]，这可能是围岩孔隙度和渗透率随深度递减的另一个原因。

5.   结论

通过热液运移和组分扩散数值模拟实验，研究硅从裂隙带向邻近围岩的扩散过程。数值模拟实验使用的流体是NaCl水溶液，其密度和粘度与温度、压力和盐度有关。本文的数值模拟实验为石英脉型钨矿床蚀变宽度垂直分带形成机制提供如下启示。

（1）温度和孔隙度在垂向上的空间分布是影响石英脉型钨矿床蚀变特征的重要变量，高温和高孔隙度有利于硅从裂隙带向邻近围岩扩散，从而形成较宽的蚀变。

（2）在围岩孔隙度均一分布的情况下，由于深部温度高于浅部，深部围岩蚀变比浅部宽，这与石英脉型钨矿床的蚀变特征相反。

（3）围岩孔隙度随深度递减会抵消温度对硅扩散的影响，可使深部围岩形成较窄的蚀变，围岩孔隙度随深度递减可能是形成石英脉型钨矿床蚀变宽度随深度减小的有效机制。相对于前人提出的岩浆热液过渡性流体假说，本文热液运移模拟为理解石英脉型钨矿床的蚀变形成过程提供了新的启示。