本项目为国家自然科学基金重点国际(地区)合作研究项目，项目名称：长江中下游成矿带典型热液矿床成矿过程三维数值模拟研究（41820104007），起止时间 2019.1-2023.12。
项目针对长江中下游成矿带热液矿床定量化的成矿作用过程及控矿因素和机制等科学问题，选择多个典型热液矿床，基于三维成矿过程数值模拟方法模拟分析了多种热液矿床成因类型（矽卡岩型、玢岩型和斑岩型）的成矿过程，以及矿床规模、矿体形态、品位分布的控制因素及其作用机理，并对热液成矿作用过程持续时间、矿床深部成矿潜力等方面开展了深入研究。
一、项目成果科普性介绍
1. 矿床学与成矿过程数值模拟
矿床学主要研究各类矿床的物质组成、成矿物质来源、成因机理及其时空分布规律，是地质科学中的主要学科之一。众多矿床中，热液矿床是最复杂、种类最多的矿床类型之一，可在不同的地质背景条件下，通过不同组成、不同来源的热液活动形成。虽然传统矿床学研究方法在成岩成矿作用、成矿机制和成矿规律等方面已有丰硕成果，但仍存在一些传统矿床学尚无法有效回答的科学问题，如成矿作用的定量化过程，矿床规模、矿体形态、品位分布及矿床定位空间的控制因素及其作用机理，成矿作用过程持续时间等。
由于热液矿床的形成过程本质上包含了一系列耦合的物理化学过程，如动力学形变、流体流动、热传导、物质运移及化学反应等（图1），因此，近年来，数值模拟方法已成为定量化研究热液矿床成矿过程的重要工具之一。
图1 热液矿床多物理化学场耦合过程
2. 典型热液矿床多场耦合数值模拟模型
通过长江中下游成矿带典型矿床的研究，项目分别建立了能够模拟矽卡岩型、玢岩型、斑岩型等类型矿床成矿过程的多场耦合数值模拟数学模型。其中的压力、传热、流体运移、成矿化学反应和物质扩散的过程控制方程如下：
(1) 传热
地质体间的热交换及整个系统的降温过程主要由方程1-5描述。
其中，T为温度，dz为简化地质几何模型的厚度，ρ为低温热液密度，Cp为热液比热容，t为计算时间，v为流体速度，q为热传导通量，q0为总流入热量，Q为总反应热，Qvd为模型与周围环境的热交换量，ρp为多孔介质（即岩体、沉积岩或煤层）密度，θp为多孔介质孔隙变化，Cp,p为含热液的多孔介质的比热容，Cpr为多孔介质导热系数，Cm为热液导热系数，keff为有效导热系数，kdisp为环境导热系数，ε为多孔介质孔隙度。
(2) 流体运移
岩体、地层和煤层中的低温热液运移过程可用达西定律进行描述（方程6）。
k为多孔介质渗透率，μ为热液粘度，g为重力常数（一般设置为9.8m/s2），ρl是流体密度。
(3) 成矿化学反应
基于质量作用定律，可用下列方程计算成矿化学反应速率。
其中，R为总反应率，vi为反应物系数，Ri为反应物反应速率，kif是正反应速率常数，ci为反应物浓度，i指代特定的反应物或生成物。
(4) 多孔介质中物质扩散
方程9、10可用于描述成矿元素在地质体间的扩散。
其中，Di是物质扩散系数，Ni为沉淀的物质的量。
3. 三维成矿过程数值模拟方法
针对热液矿床成矿过程三维数值模拟研究的实际需要，研究基于三维隐式地质建模等方法建立三维成矿过程数值模拟方法和流程。该方法能够较好的融合已有的地质资料和相关研究成果，构建用于数值计算的三维地质网格模型；基于该三维地质网格模型，能够采用有限差分或有限元的计算方法，结合定义的边界条件对上述控制方程及模型进行求解，获得瞬态或稳态的定量化模拟结果，以解析成矿过程中存在的诸多科学和矿产勘查实践问题（图2）。
图2多场耦合成矿过程数值模拟工作流程图
4. 矽卡岩型矿床成矿过程数值模拟
研究针对安庆-贵池矿集区安庆矽卡岩型铜矿床、宣城矿集区内荞麦山矽卡岩型铜矿床、麻姑山矽卡岩型铜钼矿床的成矿过程进行数值模拟。定量化描述了安庆铜矿床的双交代作用，揭示了矽卡岩化对后期矿化的作用过程，解析了蚀变分带和矿化分带的作用机理和控制因素（图3-4）；探讨了荞麦山矿床、麻姑山矿床岩体侵入体形态对成矿的控制作用，明确成矿空间分布差异的根本原因。
图3 安庆矽卡岩型铜矿床三维成矿过程数值模拟结果
图4 温度驱动下的安庆铜矿床（a）磁铁矿；（b）黄铜矿空间异位分布
研究进一步针在长江中下游成矿带安庆-贵池矿集区月山矿田和宣城矿集区麻姑山地区矽卡岩型矿床的形成过程进行了三维成矿过程数值模拟。研究发现，矽卡岩型矿床宏观上主要受到多物理场共同作用形成的扩容空间制约，是矽卡岩型矿床最为关键的定位机制和控矿条件之一；数值模拟结果不但能够用于解析矿床成矿规律和控矿机制问题，同时可为深部隐伏矿体预测能够提供重要的数据支持（图5-6）。
图5 月山矿田多场耦合三维成矿过程数值模拟与成矿预测
图6 麻姑山地区三维成矿过程数值模拟结果（VSI体应变增大率）指示深部矿化潜力
5. 玢岩型矿床成矿过程数值模拟
研究选取长江中下游成矿带宁芜矿集区内的杨庄大型铁矿床的形成过程进行了数值模拟研究。数值模拟结果显示，温度异常和接触带构造是控制杨庄铁矿床形成关键因素；杨庄铁矿床的主体形成可能只持续了1万年，化学反应作用很可能是由系统内的温度和扩散速率控制的（图7）。此外，研究综合多个IOA铁矿床的地质特征，采用多场耦合的数值模拟方法模拟了IOA铁矿床成矿系统中铁的矿化过程，探讨了IOA矿床形成的关键控制因素，约束了IOA矿床单期成矿时间（图8）。
图7 杨庄铁矿床磁铁矿化数值模拟结果
图8 不同环境模式下 IOA型矿床成矿过程（Fe3+）
（a1-a4模式I；b1-b4模式II；c1-c4模式II，阈值为1.6×10-2 mol/m3）
6. 斑岩型矿床成矿过程数值模拟
研究选取长江中下游成矿带庐枞盆地外围沙溪斑岩型铜矿床为研究对象，对其成矿过程和控矿因素进行了数值模拟研究。研究认为，围岩的物理和化学性质在斑岩成矿系统内的作用十分显著，对热液成矿系统的矿化类型、矿化强度和矿化分布具有不同的制约方式和能力。同时研究进一步认为，可以通过了解围岩孔渗系数等物理参量，预测深部隐伏斑岩型矿床的成矿深度和空间位置，相关模拟方法和结果不但具有理论意义，对于找矿勘查也具备实践应用意义（图9）。
图9 不同围岩孔深条件对沙溪斑岩型铜矿床的矿化分布和强度的影响
茶亭斑岩型铜金矿床是宣城矿集区内新发现的大型斑岩型铜-金矿床。本次工作结合矿床最新勘查资料和矿床学研究成果，通过数值模拟分析，建立了初始浓度-反应速率-矿石品位的时空关系，认为茶亭斑岩型矿床的形成过程十分短暂，可能形成于9600 – 75000年之间，与目前依据高精度同位素测年方法获得的成矿作用持续时间成果相一致。此外，基于数值模拟结果，研究对矿床-1800米以下深度的矿化潜力进行了预测，研究认为矿化可能会延续至-2400米左右，但矿化强度较浅部相比有所下降。相关预测成果可为该矿床后期的深部勘探提供支持，采用的数值模拟方法也为斑岩型矿床深部勘查提供了新的手段和方式（图10）。
图10 茶亭斑岩型铜矿床数值模拟结果
二、主要研究成果
项目通过国际合作，课题组在开放流动系统下的多物理化学场耦合三维数值模拟理论与方法和热液成矿作用过程定量模拟等方面取得创新，研究成果进一步深化了成矿作用数值模拟的理论和实践，推动和加强了长江中下游成矿带成矿理论和成矿作用研究，为深化区域成矿规律、加深对热液成矿作用方面的理解等方面提供了新的方法与依据。
依托本项目，课题组获批建立了安徽省成矿理论与找矿预测国际联合研究中心，参与本项目的国际专家进入中心开展持续性工作。另外，课题组成功获批国家留学基金委“战略性矿产资源成矿理论与找矿预测创新型人才国际合作培养项目”，将在未来选派更多优秀青年学者、研究生等前往国际知名高校进行学术交流及合作培养。
基于本项目，已发表基金标注论文21篇（其中第一标注11篇），SCI论文18篇（其中第一标注10篇），合作出版外文专著1部。