基于简单几何模型的斑岩成矿多场耦合数值模拟

Multi-Field Coupled Numerical Simulation of Porphyry Ore Formation Based on Simple Geometric Model

摘要

斑岩型矿床在全球范围内广泛分布，其成矿机制具有重要的科研价值.成矿过程数值模拟方法是研究岩浆‒热液成矿系统的重要方法，该方法能够定量、连续分析成矿作用过程，为矿化/蚀变时空分布、成矿流体迁移‒演化等问题提供解答.本研究建立了斑岩型矿床简单几何模型，并基于该模型开展了成矿过程多物理场（热‒流‒化‒质）耦合数值模拟研究.结果表明，在成矿模拟研究中使用简单几何模型具有可行性，对深部找矿靶区预测、解释大型‒超大型斑岩型矿床成因具有一定的启示意义；该方法不仅可以用来计算矿化空间分布并实现深部找矿预测，还可以通过不同的矿化分布形式来反推成矿时期的构造环境，从而更深入地研究古成矿环境等问题；此外，简单模型具有计算量少、人为影响微弱、可信度高等特点，能够在一些特定的成矿理论问题研究中发挥重要作用.

Abstract

Porphyry deposits are widely distributed in the world, and their mineralization mechanisms are of significant scientific research value. Numerical simulation is an important approach to quantitatively and continuously analyze the formation of ore deposits, revealing the temporal-spatial distribution of mineralization/alteration and migration-evolution of ore-forming fluids. This study established a simple geometric model of porphyry deposits, and conducted multi-field coupled (heat-transfer, fluid-flow, chemical reaction, diffusion) numerical simulation of its formation. Results show that using simple geometric models in metallogenic simulation research is feasible, which has certain enlightening significance for searching for deep prospecting targets and explaining the genesis of large and ultra-large porphyry deposits. This method can not only be used to calculate the spatial distribution of mineralization and achieve deep prospecting prediction, but also to infer the tectonic environment during the mineralization period through different forms of mineralization distribution, thereby enabling a deeper study of issues such as the ancient mineralization environment. In addition, simple models are characterized by low computational cost, minimal human influence, and high credibility, so they can play an important role in the study of some specific metallogenic theoretical issues.

Graphical abstract

关键词

简单几何模型 / 斑岩型矿床 / 成矿过程 / 多场耦合 / 数值模拟 / 大数据.

Key words

simple geometric model / porphyry deposit / ore formation / multi⁃field coupled / numerical simulation / big data

0 引言

斑岩型矿床是重要的铜、钼金属来源，同时还伴生有可观的金、银等贵金属，其成因与找矿勘查是矿床学领域的核心议题（Cooke et al.， 2005； Sillitoe， 2010； Weis et al.， 2012； 侯增谦等， 2020）.一般来说，斑岩型矿床规模巨大、品位低、埋藏浅，以侵入体为中心，其矿化‒蚀变系统具有显著的垂向和侧向分带性特征.系统内成矿流体以高盐度岩浆流体和后期混入的大气降水为主，金属矿物通常以浸染状和细脉状形态赋存于网状裂隙系统中.近年来，针对斑岩型矿床的研究主要聚焦于含矿岩浆的起源与演化、金属元素的超常富集机制、成矿流体出溶过程以及构造‒岩浆‒热液耦合关系等前沿问题（钟军等，2011； Mudd and Jowitt， 2018； 侯增谦等， 2020； 倪培等， 2020）.

随着数学、计算机等基础学科的发展和学科交叉融合，数学地球科学在地质学研究中逐渐发挥重要作用，地质统计学、地学大数据、计算地球科学等已被用于解决地质结构三维可视化、深部找矿预测、成矿过程定量化等方面的诸多难题（周永章等， 2018； 成秋明， 2021； 左仁广， 2025； Zuo et al.， 2025； 毛先成等， 2026）.其中，成矿过程定量化研究主要借助于多物理场耦合数值模拟方法来定量刻画成矿过程，并详细记录其中涉及的岩石形变、传热、流体迁移、成矿化学反应以及成矿物质迁移等具体物理化学过程及其内涵的复杂耦合关系，从而定量解释成矿机制、揭示深部矿化过程，为成矿机制研究和深部找矿预测提供理论支撑或直接的数据支持（Hobbs et al.， 2000； Ord et al.， 2012； Liu et al.， 2016； Zhao et al.， 2018，2022； Mao et al.， 2019； Qin et al.， 2021； Yao and Mungall， 2022）.

成矿过程数值模拟研究近年来发展迅速，随着高性能计算和多物理场耦合理论与矿床学研究相结合，当前的研究方法技术路线呈现出两个趋势：一是多过程‒多组分‒多尺度耦合的数值模拟，二是模拟与观测数据的数据同化与融合.从矿床学角度来看，当下的研究热点主要集中于成矿金属元素超常富集的关键机制（Zou et al.， 2017； Zhao et al.， 2018； Li et al.， 2019； Lai et al.， 2022； 陈伟林和肖凡， 2023； Xiong et al.， 2023； Hu et al.， 2024； Xiao et al.， 2024a，2024b）、构造‒流体成矿耦合关系（池国祥和薛春纪， 2011； Weis et al.， 2012； 李忠海等， 2015； Liu et al.， 2016； 李增华等， 2019）、“源‒运‒储”系统动态演化（Sun et al.， 2020； Lai et al.， 2022； Yao and Mungall， 2022； Leng et al.， 2025）.总的来说，这些进展推动着成矿过程数值模拟从“解释工具”向“预测工具”逐步转变，进而为成矿理论研究和深部找矿勘查提供有力支撑.

矿床尺度的成矿过程数值建模过程一般由“成矿概念模型‒简化地质模型‒数学模型”组成，简称为“三个模型”（Zou et al.， 2017）.该方法流程具有通用性与可靠性，已在众多成矿过程数值模拟研究中被成功应用（Weis et al.， 2012； Zou et al.， 2017； Zhao et al.， 2018； Lai et al.， 2022； Yao and Mungall， 2022； Hu et al.， 2024，2026）.其中，具有不同空间形态特征的简化地质模型对数值模拟结果中的流体迁移、成矿元素空间富集等结果具有强控制性，模型空间的人为区域划分和参数设定能够显著控制流体迁移轨迹和元素富集沉淀位置，对结果具有强干预性.为解决这一问题，本文尝试构建斑岩成矿系统简单地质模型，基于此开展成矿过程多场耦合数值模拟研究，并试图回答以下问题：（1）利用简单地质模型消除模型区域形态对数值模拟结果的影响后，是否还能够模拟斑岩成矿过程；（2）在基于简单地质模型的数值模拟研究中，是否能够定量描述流体迁移、元素沉淀成矿过程等关键过程.

1 研究方法

前文已述，矿床尺度成矿过程数值模拟研究一般基于“三个模型”方法流程开展，即前文所述的“成矿概念模型‒简化地质模型‒数学模型”（Zou et al.， 2017）.其中，成矿概念模型（概念建模）部分主要是对现有成矿机制的概念性、经验性描述，辅以经典矿床成矿模式图，是后续数值模拟研究的理论指导；简化地质模型（地质建模）是对研究对象的二维（剖面或水平断面）或三维（三维地质模型）简化，是开展模拟计算和结果分析的可视化载体；数学模型（数学建模）包含成矿过程中各物理化学过程的基本控制方程及其耦合关系的数学描述.随后，为模型初始、边界条件和参数、变量赋值，开展计算模拟研究与结果分析.

1.1　概念建模

斑岩成矿系统一直以来都是矿床学的研究热点，前人已对其成矿机制进行了细致的研究与总结（Sillitoe， 2010），可简述如下：在俯冲或碰撞背景下，中酸性岩浆上侵至地壳浅部；随后，从冷却结晶的岩浆中释放出富含金属和挥发分的高温流体，在岩体顶部及围岩中引发钾硅酸盐化蚀变，并沉淀出浸染状的铜钼矿体；随后，成矿流体因压力降低发生相分离，产生的酸性气体向上运移，淋滤金属并在上部形成绢英岩化和泥化蚀变带；最终，随着系统冷却，在外围形成广泛的青磐岩化带，从而构成一个以斑岩体为中心、具有典型蚀变分带和细脉浸染状矿化的完整成矿系统（图1a）.

1.2　地质建模

构建简化地质模型是本研究的核心.为了最大程度减少模型形态和区域划分对计算模拟结果的影响，基于典型斑岩‒浅成低温系统成矿模式（Sillitoe， 2010）构建了本研究的简化地质模型，模型尺寸为宽2 000 m，高2 700 m，模型空间内部无区域划分，均为岩浆固结后形成的斑岩体（图1b、1c），模型顶部深度为-1 000 m，0~-1 000 m的空间代表上覆剥蚀地层，其被模型省略，仅保留深度用以复现成矿时期的地层温度与压力.剖分后的模型含有9 412个顶点，18 474个三角形.

1.3　数学建模

本研究考虑斑岩成矿过程中的压力驱动、岩体降温、成矿流体迁移、化学反应与物质沉淀耦合过程，其数学模型可分为四部分进行描述，具体公式与符号解释如下.

1.3.1　压力驱动流体迁移

深部岩体内的成矿流体迁移过程主要受到静岩压力驱动.在已知模型空间中压力分布的情况下，可基于达西定律描述成矿流体迁移过程.此外，Kozeny⁃Carman方程可以用于描述岩石孔隙度与渗透率的动态变化（Carman， 1937）.涉及的主要控制方程如下：

P=P_{\mathrm{0}}-z\bullet G_{\mathrm{p}}，(1)

\frac{\partial }{\partial t}\left(\rho \epsilon \right)+\nabla \bullet \left(\rho \mathit{u}\right)=Q_{\mathrm{m}}，(2)

\mathit{u}=-\frac{k}{\mu }(\nabla P+\rho \mathrm{g})，(3)

k=k_{\mathrm{0}}{(\frac{\epsilon }{{\epsilon }_{\mathrm{0}}})}^{\mathrm{8}}，(4)

其中，公式（1）中，P为压力因变量（单位：Pa），P_{\mathrm{0}}是地表压力（单位：Pa），z是深度（单位：m），G_{\mathrm{p}}为压力梯度（单位：Pa/m）；公式（2）中，t为时间（单位：s），\rho为密度（单位：kg/m³），\epsilon为孔隙度因变量，\mathit{u}为流体速度矢量场（单位：m/s），Q_{\mathrm{m}}为质量源（单位：kg/（m³‧s））；公式（3）与（4）中，k为渗透率（单位：m²），\mu为成矿流体黏度（单位：Pa‧s），\mathrm{g}为重力常数（取值为9.8 m/s²），k_{\mathrm{0}}为岩石初始渗透率（单位：m²），{\epsilon }_{\mathrm{0}}为岩石初始孔隙度.

1.3.2　岩体降温

岩石内部的热变化方式主要为热传递，在岩浆‒热液系统成矿过程中具有核心地位，驱动着降温活动并控制岩体内部成矿化学反应的进行（Fournier， 1999）.根据模型的深度、地表温度以及地温梯度（胡宝群等， 2003），笔者使用以下方程计算模型的初始温度：

T=T_{\mathrm{0}}-z·G_T,(5)

其中，T为模型各处的温度因变量，T_{\mathrm{0}}为地表温度（设置为20 °C，压力为一个标准大气压，取值为101.325 kPa），G_{\mathrm{T}}为温度梯度（取值为0.026 °C/m）.此外，模型中的热过程（能量交换、温度变化等）用以下方程定量描述：

d_{\mathrm{z}}(\rho C_{\mathrm{p}}{)}_{\mathrm{e}\mathrm{f}\mathrm{f}}\frac{\partial T}{\partial t}+d_{\mathrm{z}}{\rho }_{\mathrm{f}}{C}_{\mathrm{p},\mathrm{f}}\mathit{u}\bullet \nabla T+\nabla \bullet \mathit{q}=d_{\mathrm{z}}Q+q_{\mathrm{0}}+d_{\mathrm{z}}{Q}_{\mathrm{p}}+d_{\mathrm{z}}{Q}_{\mathrm{v}\mathrm{d}}+d_{\mathrm{z}}{Q}_{\mathrm{g}\mathrm{e}\mathrm{o}}，(6)

\mathit{q}=-d_z{\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}\nabla T，(7)

(\rho C_{\mathrm{p}}{)}_{\mathrm{e}\mathrm{f}\mathrm{f}}={\epsilon }_{\mathrm{s}}{\rho }_{\mathrm{s}}{C}_{\mathrm{p},\mathrm{s}}+{\epsilon }_{\mathrm{f}}{\rho }_{\mathrm{f}}{C}_{\mathrm{p},\mathrm{f}}+{\epsilon }_{\mathrm{i}\mathrm{m}\mathrm{t}}{\rho }_{\mathrm{i}\mathrm{m}\mathrm{t}}{C}_{\mathrm{p},\mathrm{i}\mathrm{m}\mathrm{t}}，(8)

{\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}={\epsilon }_{\mathrm{s}}{\lambda }_{\mathrm{s}}+{\epsilon }_{\mathrm{p}}{\lambda }_{\mathrm{f}}+{\epsilon }_{\mathrm{i}\mathrm{m}\mathrm{t}}{\lambda }_{\mathrm{i}\mathrm{m}\mathrm{t}}+{\lambda }_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p}}，(9)

其中，公式（6）中，d_{\mathrm{z}}为模型厚度（单位：m），C_{\mathrm{p}}是含有成矿流体的地质体综合比热容（单位：J/（kg‧K）），(\rho C_{\mathrm{p}}{)}_{\mathrm{e}\mathrm{f}\mathrm{f}}为有效恒压热容（单位：J/（kg‧K），下标eff为effective，即有效部分），{\rho }_{\mathrm{f}}为流体密度（单位：g/m³，下标f为fluid，即流体），C_{\mathrm{p},\mathrm{f}}为流体比热容（单位：J/（kg‧K）），\mathit{u}为流体速度矢量场（单位：m/s），\mathit{q}为热传导通量矢量（单位：W/m²），Q为全局热源/广义内热源（单位：W/m³），q_{\mathrm{0}}为平面/线热源（单位：W/m²），Q_{\mathrm{p}}为压力驱动流体迁移产生的黏性耗散热源（单位：W/m³），Q_{\mathrm{v}\mathrm{d}}为体积形变对模型系统做功（单位：W/m³），Q_{\mathrm{g}\mathrm{e}\mathrm{o}}为背景的、恒定的地热源（单位：W/m³）；公式（7）中，{\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}为有效导热率；公式（8）与（9）中，\epsilon为体积分数（对于岩石来说，体积分数{\epsilon }_{\mathrm{s}}=\mathrm{1}-\epsilon，即岩石中固体部分所占的比率，而对于岩石孔隙中的流体来说，体积分数{\epsilon }_{\mathrm{f}}=\epsilon，即岩石孔隙度本身），{\rho }_{\mathrm{s}}为岩石固体密度（单位：kg/m³），C_{\mathrm{p},\mathrm{s}}为岩石固体比热容（单位：J/（kg‧K）），C_{\mathrm{p},\mathrm{f}}为流体比热容（单位：J/（kg‧K））；下标为imt的{\epsilon }_{\mathrm{i}\mathrm{m}\mathrm{t}}、{\rho }_{\mathrm{i}\mathrm{m}\mathrm{t}}、C_{\mathrm{p},\mathrm{i}\mathrm{m}\mathrm{t}}分别为中间相/夹杂物的体积分数、密度以及比热容，其并不属于多孔介质的固体框架，也不属于流体填充物，但{\epsilon }_{\mathrm{i}\mathrm{m}\mathrm{t}}\mathrm{、}{\rho }_{\mathrm{i}\mathrm{m}\mathrm{t}}\mathrm{、}{C}_{\mathrm{p},\mathrm{i}\mathrm{m}\mathrm{t}}项的存在使模型能够适用于相对复杂的情形.

1.3.3　化学反应与元素沉淀

本研究以铜元素为例（斑岩型矿床中最主要的成矿金属元素），使用简化的化学反应方程描述铜元素沉淀成矿的过程，以斑岩阶段黄铜矿生成的温度区间作为成矿温度区间，即550~400 °C（单思齐等， 2023），其反应式如下：

\mathrm{C}{\mathrm{u}}^{\mathrm{2}+}\stackrel{R}{\to }\mathrm{C}\mathrm{u}\downarrow，(10)

其中，R表示该过程经历了化学反应.与此同时，质量作用定律能够用于计算反应物/生成物的增加/减少过程中的动态化学反应速率（Grmela， 2012）：

R={\sum }_i{I}_i{r}_i，(11)

r_i=k_i^j{\prod }_{i∊\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{c}\mathrm{t}}{c}_i^{-I_i}，(12)

其中，R为全局化学反应速率（单位：mol/（m³‧s）），I_i为各反应物/生成物的计量系数，r_i为物质i的反应速率（单位：mol/（m³‧s）），k_i^j为正反应速率常数，j为正反应次数，c_i是流体中物质i的浓度（单位：mol/m³）.

由于直接用生成物的累计量描述物质沉淀并不能完全移除物质扩散对沉淀累计的影响，因此，物质沉淀无法使用反应物/生成物的化学反应方程式进行直接描述.本研究基于值传递方法，使用分布式域常微分方程在模型各处重新计算元素沉淀生成的积分值，方程如下：

f=e_a\frac{{\partial }^{\mathrm{2}}{c}_i}{\partial t^{\mathrm{2}}}+d_a\frac{\partial c_i}{\partial t}，(13)

其中，f为微分方程源项，此处为r_{\mathrm{C}\mathrm{u}}，即铜元素的反应速率（单位：mol/（m³‧s）），e_a为质量系数，理论上为0，但一般设置为一个极小数（如1×10^-15）来加速计算收敛，d_a为阻尼系数，设置为1.模拟过程中，对源项f进行域积分，即可计算获得模型域空间内各处积累的且不受扩散影响的生成物的总量.

1.3.4　物质扩散

物质扩散是在小空间尺度上影响矿化空间分布的重要过程之一.下列方程主要用来模拟模型空间中的\mathrm{C}{\mathrm{u}}^{\mathrm{2}+}物质扩散过程：

\frac{\partial c_i}{\partial t}+\nabla \bullet \left(-D_i\nabla c_i\right)+\mathrm{\nu }\bullet \nabla c_i=r_i，(14)

N_i=-D_i\nabla c_i+\nu c_i，(15)

其中，c_i是反应物/生成物i的物质浓度（单位：mol/m³），D_i是物质i的扩散系数（单位：m²/s），\mathrm{\nu }为模型空间内的流体速度（单位：m/s），r_i是对应的物质i的化学反应速率（单位：mol/（m³‧s）），N_i是这一过程中沉淀的物质i的总量（单位：mol）.

1.4　模拟设定

丰富的前人研究为成矿模拟积累了坚实的理论基础与大量可用数据（Weis et al.， 2012； Hu et al.， 2020； 单思齐等， 2023）.根据前人绘制的矿物生成顺序表，设置成矿化学反应温度区间为550~400 °C，整个岩体的初始温度被设置为 600 °C，以便于整个模型的降温过程能够完全涵盖成矿化学反应温度区间，以增加合理性；模型的初始温度和压力根据模型深度以及常用的地温梯度/地压梯度计算，模型左右两侧的压力均设置为向内的8×10⁷ Pa，主要目的是通过施加压力构建等价的构造应力环境，同时实现流体汇聚；考虑到质量作用定律的限制，将成矿流体中携带的\mathrm{C}{\mathrm{u}}^{\mathrm{2}+}的浓度设置为0.05 mol/m³；根据单一降温过程数值模拟结果中的降温持续时间，将数值模拟时限设置为60 000年，含矿流体出现的时间段为1 000~10 000年，其余时间从底部向上迁移的流体为非含矿流体，以模拟单期成矿过程；时间步长设置为10年，并生成瞬态结果；岩体是模型空间内部仅有的唯一地质单元，其岩石物性参数见表1.

2 结果

2.1　降温

本研究中模型为整块岩体中的一部分，可视为均质模型，因此，降温过程表现较为平稳，整个模型范围（岩体）由初始的600 °C逐步缓慢降温，直至温度降低至成矿化学反应温度区间以下，即550~400 °C.期间，适宜成矿化学反应发生的域内温度区间平稳下移，直至最终离开模型底部（图2）.

2.2　流体携成矿金属元素迁移

由于受到模型左右两侧的向内挤压，成矿流体的迁移方向较为复杂（图3a），模型中部形成了向上的流体迁移高速区（图3b）.受单期含矿流体 （1 000~10 000万年）的影响，模型底部自10 000年后，\mathrm{C}{\mathrm{u}}^{\mathrm{2}+}的浓度逐渐降低，而模型内部残存的含矿流体逐步向上迁移，从之前形成的流体迁移高速区流出，部分成矿金属元素在成矿温度区间（550~400 °C）发生反应并沉淀成矿；自30 000年后，模型空间内仅存在少量\mathrm{C}{\mathrm{u}}^{\mathrm{2}+}，并在后续的 30 000年间逐步随流体迁移离开模型空间（图4）.

2.3　成矿化学反应

受到降温过程影响，成矿过程化学反应仅在550~400 °C发生，该温度区间与流体迁移中的高速区在模型中段相遇，反应结果可用公式（12）所计算的域积分结果表达（图5）.显然，降温过程中550~400 °C温度区间与成矿金属元素随流体向上迁移的交集的动态范围即是最适宜成矿元素卸载沉淀成矿范围，发生沉淀的Cu元素最大浓度达到了5×10³ mol/m³，相对于成矿金属元素迁移过程中的最大浓度5×10^-2 mol/m³，其富集程度达10⁵倍.随后，改变两侧压力值进行一组对比实验，左侧向内压力值保持为8×10⁷ Pa，右侧向内压力改为7.5× 10⁷ Pa，双侧压力值略有差异，目的是打破轴对称，对模拟结果增加一定程度的非规则性，以贴近现实情况中地质环境的随机现象.计算结果显示，矿化的空间分布发生了一定程度的侧向偏移（图6）.

3 讨论

本研究建立了斑岩型矿床简单模型，并基于此开展了多物理场耦合成矿过程数值模拟研究.由于使用了简单几何模型，模拟结果受到了尽可能少的人为域空间划分影响，但在合理的条件设置下，模拟结果仍然表现出了斑岩型矿床成矿过程的显著特征，在不同的侧向压力下，模型内部的铜矿化分布分别展现出了部分斑岩型矿床的矿化空间分布特征.上述研究结果表明，在成矿模拟研究中使用简单几何模型是可行性，相较于使用复杂划分的几何模型，简单几何模型能在减小人为干预因素之后，表现出更具可信度的模拟结果.

从矿化形成的角度来看，本研究基于简单几何模型，耦合了斑岩系统成矿过程中的降温、流体迁移、成矿化学反应、物质扩散等几个主要物理、化学过程，这些条件均对矿床形成过程具有一定的约束与限制，只有在同时满足所有适宜成矿的条件时，才能够在模型空间的特定位置形成矿化累积，这对于理解斑岩型矿床的成矿机制具有重要意义，同时也解释了为何世界上大型、超大型矿床较为稀有，一旦发现就会成为矿床学研究热点（Meinert et al.， 2005； Sillitoe， 2010）.

从流体迁移的角度来看，本研究中流体携带成矿金属元素的向上迁移过程受到了侧向压力的极大影响，其迁移轨迹受压力、浮力的共同作用.与此同时，两侧压力（或其他等价构造条件）是造成斑岩成矿过程中流体汇聚的必要条件，而侧向压力值的差异则是形成不同形态矿化分布的主要控制条件.如两侧压力相近或相等，则会形成矿化垂向分布的斑岩型矿床，矿化分布从下向上近似垂直，如位于美国犹他州的Bingham、内华达州的Yerington矿床、中国安徽省的茶亭矿床均为此类（Weis et al.， 2012； Hu et al.， 2020； 图6）；如果两侧压力相差较大，则会形成另一类典型的斑岩型矿床，矿化从下向上有一定的偏移，我国福建紫金山斑岩‒浅成低温系统属于此类（王少怀等， 2009； 钟军等， 2011； 图7）.因此，成矿模拟方法不仅可以用来计算矿化空间分布并实现找矿预测，还可以通过不同的矿化展布形式来反推成矿时期的构造环境，从而更深入地研究成矿古环境等问题.

从简单几何模型的角度来看，本文尽可能减少了模型内部空间的几何划分，从而最大程度减少了人为干预因素.相对于有着复杂划分的几何模型，简单几何模型有以下优势：（1）相对计算量较 少——在常规的二维/三维成矿过程数值模拟研究中，模型空间往往被划分为数个区域，这些区域在其边缘往往有较为细碎的模型剖分结果，这在建模过程、计算过程中有概率会导致不稳定性，造成计算量大、模型不收敛等问题，影响研究进程；（2）更凸显成矿条件的重要性——相对于有着复杂划分的几何模型来说，简单几何模型只有一套岩石物性参数，其内部传热、流体迁移、成矿化学反应等过程均需要综合设置各项条件后才能表现得接近真实情况，而对这些成矿条件（如侧向压力、流体迁移轨迹、P⁃T⁃x曲线等）进行细致调节，也是逐步反推古成矿环境的过程.相比之下，如果使用复杂划分的几何模型，则需要分别为不同的域空间设置不同的岩石物性参数，降温、流体迁移、成矿化学反应等过程就会极大地受到人为因素干预（Hu et al.， 2024），导致模型的可信度相对较低，难以解释复杂的成矿作用.例如，前期针对安徽宣城茶亭斑岩型铜矿床开展的成矿模拟研究使用了复杂剖分的几何模型，其矿化形成过程及其空间展布受到了较强的主观控制影响（Hu et al.， 2020），而针对安徽月山朱冲富钴矽卡岩型铁矿床开展的成矿模拟研究则使用了较为简单的地质模型，其降温过程使铁矿体中的钴出现了规律的“条带型”富集（Hu et al.， 2026），对于定量剖析其成矿作用具有启示意义（图7）.

然而，本研究也存在一些不足之处.首先，是依旧未能更精细地描述成矿化学反应过程，简化的化学反应方程式能够减少计算量，在研究中快速发挥作用（Grmela， 2012），但与真实的成矿化学反应组相比，依旧相去甚远，如要研究精细的成矿过程，仍旧需要建立精细地、耦合的成矿化学反应方程组，而非简化方程.其次，基于简单几何模型的成矿模拟只能解决部分矿床的部分问题，如果研究对象发生改变，则难以适应.举例来说，在斑岩型矿床成矿数值模拟研究中，可以取岩体的一部分建立一个形态规则的简单几何模型，这是因为斑岩型矿床具有岩体范围大且相对均一的形态特征，如果研究对象是矽卡岩型矿床，其矿化主要形成在岩体与地层的接触带，且空间分布范围狭小，而岩体和地层这两类地质单元的岩石物性参数相差巨大，建模时必须进行几何划分，不能将整个模型划归为一个整体.此外，对于具体的矿床学问题而言，基于简单几何模型能够对矿化展布特征、古成矿环境、主要控矿因素等问题作出解答，但对于成矿持续时间等问题则难以给出答案（Xiao et al.， 2024b）.除此之外，经典的成矿过程数值模拟工作流程虽然能够解决一部分矿床学问题，但仍面临参数优化、模拟结果与地质观测事实吻合度差等问题，如能结合机器学习、数据同化等方法，将有可能解决上述问题，从而获得更精确、更可信的模拟结果.例如，基于数据同化理论，通过迭代优化不断校正模拟参数与边界条件，进行“数值模拟‒参数优化‒数据同化”的计算循环，能够实现模拟参数优选并使得模拟结果尽可能逼近地质事实，从而生成更精准的成矿预测数据（如矿化空间富集区、成矿流体高通量区等），为成矿预测研究提供有利信息.总而言之，未来随着数值模拟理论、方法、计算机软硬件、分析测试方法的进步以及机器学习、数据同化在数值模拟研究中的应用，这些科学问题都将逐步得到解决并推动成矿模拟领域进一步发展.

4 结论

本研究基于简单模型和多物理场耦合方法，对斑岩型矿床成矿过程开展了数值模拟研究，得到了以下结论：

（1）模拟结果验证了在成矿模拟研究中使用简单几何模型的可行性，相较于复杂划分几何模型，简单模型具有计算量少、人为影响微弱、可信度高等特点，能够在部分成矿理论问题研究中发挥作用.

（2）成矿模拟方法不仅可以用来计算矿化空间分布并实现深部找矿预测，还可以通过不同的矿化展布形式来反推成矿时期的构造环境，从而更深入地研究古成矿环境等问题.

（3）当前，成矿过程数值模拟研究仍存在数学模型理想化/过于简化、泛用性低、难以解决部分矿床学难题（如成矿持续时间）等一系列问题，但未来随着数值模拟理论、方法、计算机软硬件、分析测试方法的进步以及机器学习、数据同化在数值模拟研究中的应用，这些科学问题都将逐步得到解决并推动成矿模拟领域进一步发展.

参考文献

[1] Carman, P. C., 1937. Fluid Flow Through Granular Beds. Transactions of the Institution of Chemical Engineers, 15: 150-166. https://doi.org/10.1016/S0263⁃8762(97)80003⁃2

[2] Chen, W. L., Xiao, F., 2023. Advances in Numerical Modeling of Metallogenic Dynamics: A Review of Theories, Methods and Technologies. Bulletin of Geological Science and Technology, 42(3): 234-249 (in Chinese with English abstract).

[3] Cheng, Q. M., 2021. What Are Mathematical Geosciences and Its Frontiers? Earth Science Frontiers, 28(3): 6-25 (in Chinese with English abstract).

[4] Chi, G. X., Xue, C. J., 2011. Principles, Methods and Applications of Hydrodynamic Studies of Mineralization. Earth Science Frontiers, 18(5): 1-18 (in Chinese with English abstract).

[5] Cooke, D. R., Hollings, P., Walshe, J. L., 2005. Giant Porphyry Deposits: Characteristics, Distribution, and Tectonic Controls. Economic Geology, 100(5): 801-818. https://doi.org/10.2113/gsecongeo.100.5.801

[6] Fournier, R. O., 1999. Hydrothermal Processes Related to Movement of Fluid from Plastic into Brittle Rock in the Magmatic⁃Epithermal Environment. Economic Geology, 94(8): 1193-1211. https://doi.org/10.2113/gsecongeo.94.8.1193

[7] Grmela, M., 2012. Fluctuations in Extended Mass⁃Action⁃Law Dynamics. Physica D: Nonlinear Phenomena, 241(10): 976-986. https://doi.org/10.1016/j.physd.2012.02.008

[8] Hobbs, B. E., Zhang, Y., Ord, A., et al., 2000. Application of Coupled Deformation, Fluid Flow, Thermal and Chemical Modelling to Predictive Mineral Exploration. Journal of Geochemical Exploration, 69: 505-509. https://doi.org/10.1016/S0375⁃6742(00)00099⁃6

[9] Hou, Z. Q., Yang, Z. M., Wang, R., et al., 2020. Further Discussion on Porphyry Cu⁃Mo⁃Au Deposit Formation in Chinese Mainland. Earth Science Frontiers, 27(2): 20-44 (in Chinese with English abstract).

[10] Hu, B. Q., Wang, F. Z., Sun, Z. X., et al., 2003. The Pressure Gradient in the Lithosphere. Earth Science Frontiers, 10(3): 129-133 (in Chinese with English abstract).

[11] Hu, X. Y., Chen, R. Y., Liu, G. X., et al., 2024. Numerical Simulation⁃Based Targeting of the Au Mineralization within the Hongbu Altered Rock Type Gold Deposit, Eastern Shandong Peninsula, China. Journal of Geochemical Exploration, 265: 107557. https://doi.org/10.1016/j.gexplo.2024.107557

[12] Hu, X. Y., Li, X. H., Yuan, F., et al., 2020. Numerical Modeling of Ore⁃Forming Processes within the Chating Cu⁃Au Porphyry⁃Type Deposit, China: Implications for the Longevity of Hydrothermal Systems and Potential Uses in Mineral Exploration. Ore Geology Reviews, 116: 103230. https://doi.org/10.1016/j.oregeorev.2019.103230

[13] Hu, X. Y., Liang, X., Wang, F. Y., et al., 2026. Cooling Front Controls on Co Mineralization within Skarn⁃Type Co⁃Rich Fe Deposits: A Numerical Modeling Case Study of the Zhuchong Deposit, Anhui Province, China. Natural Resources Research, 35(1): 17-31. https://doi.org/10.1007/s11053⁃025⁃10573⁃z

[14] Lai, F., Zou, S. H., Xu, D. R., 2022. Silver Complexation in Chlorine⁃ and Sulfur⁃Rich Hydrothermal Fluids: Insight from Ab Initio Molecular Dynamics Simulations. Chemical Geology, 589: 120684. https://doi.org/10.1016/j.chemgeo.2021.120684

[15] Leng, S., Xia, Q. L., Li, T. F., et al., 2025. Geodynamic Simulation of the Pulang Porphyry Deposit in Southwest China: Implications for Ore Genesis and Exploration. Journal of Earth Science, 36(6): 2513-2525. https://doi.org/10.1007/s12583⁃022⁃1727⁃9

[16] Li, X. H., Yuan, F., Zhang, M. M., et al., 2019. 3D Computational Simulation⁃Based Mineral Prospectivity Modeling for Exploration for Concealed Fe⁃Cu Skarn⁃Type Mineralization within the Yueshan Orefield, Anqing District, Anhui Province, China. Ore Geology Reviews, 105: 1-17. https://doi.org/10.1016/j.oregeorev.2018.12.003

[17] Li, Z. H., Chi, G. X., Deng, T., et al., 2019. Controls of Reactivated Faults on the Unconformity⁃Related Uranium Deposits in the Athabasca Basin, Canada. Geotectonica et Metallogenia, 43(3): 518-527 (in Chinese with English abstract).

[18] Li, Z. H., Liu, M. Q., Gerya, T., 2015. Material Transportation and Fluid⁃Melt Activity in the Subduction Channel: Numerical Modeling. Scientia Sinica (Terrae), 45(7): 881-899 (in Chinese with English abstract).

[19] Liu, L. M., Li, J. F., Zhou, R. C., et al., 2016. 3D Modeling of the Porphyry⁃Related Dawangding Gold Deposit in South China: Implications for Ore Genesis and Resources Evaluation. Journal of Geochemical Exploration, 164: 164-185. https://doi.org/10.1016/j.gexplo.2015.11.002

[20] Mao, X. C., Duan X. M., Deng, H., et al., 2026. Intelligent 3D Prediction of Deep Mineral Resources: Theory, Methods, and Challenges. Earth Science, 51(3): 793-815 (in Chinese with English abstract).

[21] Mao, X. C., Ren, J., Liu, Z. K., et al., 2019. Three⁃ Dimensional Prospectivity Modeling of the Jiaojia⁃Type Gold Deposit, Jiaodong Peninsula, Eastern China: A Case Study of the Dayingezhuang Deposit. Journal of Geochemical Exploration, 203: 27-44. https://doi.org/10.1016/j.gexplo.2019.04.002

[22] Meinert, L. D., Dipple, G. M., Nicolescu, S., 2005. World Skarn Deposits. Economic Geology, (One Hundredth Anniversary Volume): 299-336. https://doi.org/10.5382/AV100.11

[23] Mudd, G. M., Jowitt, S. M., 2018. Growing Global Copper Resources, Reserves and Production: Discovery Is not the Only Control on Supply. Economic Geology, 113(6): 1235-1267. https://doi.org/10.5382/econgeo.2018.4590

[24] Ni, P., Chi, Z., Pan, J. Y., 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).

[25] Ord, A., Hobbs, B. E., Lester, D. R., 2012. The Mechanics of Hydrothermal Systems: I. Ore Systems as Chemical Reactors. Ore Geology Reviews, 49: 1-44. https://doi.org/10.1016/j.oregeorev.2012.08.003

[26] Qin, Y. Z., Liu, L. M., Wu, W. C., 2021. Machine Learning⁃Based 3D Modeling of Mineral Prospectivity Mapping in the Anqing Orefield, Eastern China. Natural Resources Research, 30(5): 3099-3120. https://doi.org/10.1007/s11053⁃021⁃09893⁃7

[27] Shan, S. Q., Xie, G. Q., Liu, W. Y., et al., 2023. Hydrothermal Alteration and Sulfide Zoning Characteristics of the Čukaru Peki Super⁃Large Porphyry⁃Epithermal Copper⁃Gold Deposit in Western Tethys, Serbia and Its Implication for Exploration. Geotectonica et Metallogenia, 47(5): 1085-1109 (in Chinese with English abstract).

[28] Sillitoe, R. H., 2010. Porphyry Copper Systems. Economic Geology, 105(1): 3-41. https://doi.org/10.2113/gsecongeo.105.1.3

[29] Sun, T., Li, H., Wu, K. X., et al., 2020. Data⁃Driven Predictive Modelling of Mineral Prospectivity Using Machine Learning and Deep Learning Methods: A Case Study from Southern Jiangxi Province, China. Minerals, 10(2): 102. https://doi.org/10.3390/min10020102

[30] Wang, S. H., Pei, R. F., Zeng, X. H., et al., 2009. Metallogenic Series and Model of the Zijinshan Mining Field. Acta Geologica Sinica, 83(2): 145-157 (in Chinese with English abstract).

[31] Weis, P., Driesner, T., Heinrich, C. A., 2012. Porphyry⁃Copper Ore Shells Form at Stable Pressure⁃Temperature Fronts within Dynamic Fluid Plumes. Science, 338(6114): 1613-1616. https://doi.org/10.1126/science.1225009

[32] Xiao, F., Chen, X. Y., Cheng, Q. M., 2024a. Combining Numerical Modeling and Machine Learning to Predict Mineral Prospectivity: A Case Study from the Fankou Pb⁃Zn Deposit, Southern China. Applied Geochemistry, 160: 105857. https://doi.org/10.1016/j.apgeochem.2023.105857

[33] Xiao, F., Wang, K. Q., Cheng, Q. M., 2024b. Porphyry Magma Cooling and Crystallization Control of Mineralization: Insights from the Dynamic Numerical Modeling. Ore Geology Reviews, 166: 105956. https://doi.org/10.1016/j.oregeorev.2024.105956

[34] Xiong, Y. H., Zuo, R. G., Miller, S. A., 2023. The Behavior of Hydrothermal Mineralization with Spatial Variations of Fluid Pressure. Journal of Geophysical Research: Solid Earth, 128(2): e2022JB025255. https://doi.org/10.1029/2022JB025255

[35] Yao, Z. S., Mungall, J. E., 2022. Magnetite Layer Formation in the Bushveld Complex of South Africa. Nature Communications, 13: 416. https://doi.org/10.1038/s41467⁃022⁃28000⁃9

[36] Zhao, C. B., Hobbs, B. E., Ord, A., 2018. Modeling of Mountain Topography Effects on Hydrothermal Pb⁃Zn Mineralization Patterns: Generic Model Approach. Journal of Geochemical Exploration, 190: 400-410. https://doi.org/10.1016/j.gexplo.2018.04.004

[37] Zhao, C. B., Hobbs, B. E., Ord, A., 2022. Investigating Porosity Structure Effects on Hydrothermal Mineralization Patterns within Permeable Horizontal Layers of Fluid⁃Saturated Rocks: Semi⁃Analytical Approach through Generic Models. Ore Geology Reviews, 150: 105116. https://doi.org/10.1016/j.oregeorev.2022.105116

[38] Zhong, J., Chen, Y. J., Chen, J., et al., 2011. Fluid Inclusion Study of the Luoboling Porphyry Cu⁃Mo Deposit in the Zijinshan Ore Field, Fujian Province. Acta Petrologica Sinica, 27(5): 1410-1424 (in Chinese with English abstract).

[39] Zhou, Y. Z., Chen, S., Zhang, Q., et al., 2018. Advances and Prospects of Big Data and Mathematical Geoscience. Acta Petrologica Sinica, 34(2): 255-263 (in Chinese with English abstract).

[40] Zou, Y. H., Liu, Y., Dai, T. G., et al., 2017. Finite Difference Modeling of Metallogenic Processes in the Hutouya Pb⁃Zn Deposit, Qinghai, China: Implications for Hydrothermal Mineralization. Ore Geology Reviews, 91: 463-476. https://doi.org/10.1016/j.oregeorev.2017.09.008

[41] Zuo, R. G., Yang, F. F., Cheng, Q. M., et al., 2025. A Novel Data⁃Knowledge Dual⁃Driven Model Coupling Artificial Intelligence with a Mineral Systems Approach for Mineral Prospectivity Mapping. Geology, 53(3): 284-288. https://doi.org/10.1130/g52970.1

[42] Zuo, R. G., 2025. Key Technology for Intelligent Mineral Prospectivity Mapping: Challenges and Solutions. Scientia Sinica (Terrae), 55(9): 3104-3119 (in Chinese with English abstract).

基金资助

国家自然科学基金项目(42202327)

国家自然科学基金项目(42230802)

科技部深地重大专项课题(2024ZD1002100)

科技部深地重大专项课题(2024ZD1001404)

安徽省战略性矿产资源深部探测与评价利用安徽省重点实验室开放基金项目(PA2025GDSK0086)