关键词：反应输运模拟、热力学平衡、PHREEQC、TOUGHREACT、CrunchFlow、砂岩型铀矿、地下水化学演化、氧化还原前锋、蒙特卡洛不确定性分析、同位素混合模型、吉布斯自由能最小化、质量作用定律、工作站配置
技术标签：反应-输运模拟、热力学平衡计算、蒙特卡洛不确定性分析、同位素混合模型、PHREEQC、TOUGHREACT、CrunchFlow、吉布斯自由能最小化、质量作用定律、氧化还原反应、地下水流动、溶质运移、数值求解、参数敏感性分析
一、从”铀随水去”到”水去铀留”：一场地下水里的氧化还原接力
1983年1月，北京大学地质系举办了一场看似不起眼的讲习班——”稳定同位素地球化学基础理论讲习班”。主讲人郑淑慧副教授，听众来自全国三十多个单位的四十多名科技人员。这是中国第一次系统讲授氢氧同位素分馏、矿物同位素地质温度计、岩浆分异过程的同位素特征。
没人想到，这场讲习班播下的种子，会在四十年后长成一片森林。今天的中国地球化学家，用同位素示踪地下水循环、用反应-输运模拟预测污染物迁移、用蒙特卡洛量化成矿不确定性——所有这些工具的底层逻辑，都可以追溯到1983年北大那间教室里的板书。
而其中最”硬核”的应用之一，是砂岩型铀矿的成矿预测。
砂岩型铀矿是中国铀资源最重要的类型，占已探明储量的60%以上。它的成矿机制听起来像一句顺口溜：“铀随水去，水去铀留”。
前半句”铀随水去”——含氧地下水流经蚀源区，把岩石里的四价铀（UO₂，不溶于水）氧化成六价铀（UO₂²⁺，易溶于水），然后随着地下水向盆地中心迁移。后半句”水去铀留”——当含铀水流遇到还原环境（有机质、硫化氢、油气流体等），六价铀被还原回四价铀，以沥青铀矿（UO₂）的形式沉淀下来，富集成矿。
这个氧化-还原界面，就是砂岩型铀矿的地球化学前锋线，也叫卷头（roll front）。1956年美国怀俄明盆地首次发现这种”卷状”矿体时，地质学家是靠钻孔取样和手绘剖面来推断前锋线位置的。今天，我们用三维反应-输运模拟在计算机里”重演”整个成矿过程——从百万年前的地下水补给，到氧化还原前锋的推进，再到铀矿体的最终定位。
但这里的算力需求，远比大多数人想象的复杂。
二、反应-输运模拟：地下水化学演化的数值重现
要理解砂岩型铀矿的形成，必须同时解决三个物理化学过程：地下水流动、溶质运移、化学反应。这三个过程互相耦合：水流驱动溶质迁移，化学反应改变溶液成分和矿物相，矿物溶解/沉淀又改变孔隙度和渗透系数，反过来影响水流。
这就是**反应-输运模拟（Reactive Transport Modeling, RTM）**的核心挑战。
2.1 控制方程的本质
反应-输运模拟的数学骨架由三个方程组成：
地下水流动方程（Richards方程或达西定律）：描述多孔介质中水的运动，受渗透系数、水力梯度、孔隙度控制。对于非饱和带，还要考虑毛细压力和相对渗透率的非线性关系。
溶质运移方程（对流-弥散方程）：描述溶解组分的迁移，包括对流（随水流运动）、水动力弥散（机械弥散+分子扩散）、吸附/解吸、放射性衰变。对于铀矿模拟，还要考虑六价铀和四价铀的氧化还原转化。
化学反应方程：描述水-岩-气三相之间的化学反应，包括矿物溶解/沉淀、离子交换、表面络合、氧化还原反应、气体溶解/逸出。这些反应可以是瞬时平衡（热力学控制）或动力学控制（反应速率限制）。
这三个方程的耦合方式决定了数值求解的复杂度。最常见的策略是算子分裂（Operator Splitting）：在每个时间步内，先求解流动方程得到流速场，再求解运移方程得到组分分布，最后求解化学反应方程更新矿物饱和度和溶液成分。分裂的精度取决于时间步长——步长太大，化学反应和运移的耦合误差会累积；步长太小，计算时间会爆炸。
2.2 主流软件对比
目前地球化学模拟领域有三款主流代码，各有侧重：
PHREEQC（美国地质调查局USGS开发）：地球化学模拟的”瑞士军刀”。核心功能是批次反应计算（化学平衡、反应路径、反向模拟）和一维运移-反应耦合。PHREEQC的化学引擎极其强大，支持离子缔合模型（如Debye-Hückel、Pitzer）、多种表面络合模型（扩散双层、CD-MUSIC）、动力学反应模块。但它的运移模块仅限一维，三维场景需要耦合外部流动模型（如MODFLOW-PHAST、OGS-PHREEQC）。
TOUGHREACT（劳伦斯伯克利国家实验室开发）：基于TOUGH2多相流模拟器，专攻非等温多相流体-反应耦合。它用有限体积法处理复杂几何和不规则网格，支持双孔隙/双渗透率模型（裂隙-基质交换）、多区域多相互作用连续体（MINC）。TOUGHREACT的化学反应模块相对简单，主要处理平衡反应和一级动力学，但它在处理地热系统、CO₂地质封存、核废料处置等多相流+热+化耦合问题时几乎无可替代。
CrunchFlow（ also 劳伦斯伯克利国家实验室）：专攻多组分反应-输运，采用全耦合（全局隐式）求解策略而非算子分裂，在强耦合反应（如矿物沉淀堵塞孔隙）时精度更高。CrunchFlow支持自适应网格细化（AMR），在反应前锋附近自动加密网格，节省计算量。
对于砂岩型铀矿成矿模拟，选择哪款软件取决于研究目标：
- PHREEQC适合一维柱实验模拟、反向地球化学建模（从观测水化学推断反应路径）、快速参数扫描
- TOUGHREACT适合盆地尺度三维模拟，特别是涉及非饱和带、地热梯度、多相流（油气-水-气）的复杂场景
- CrunchFlow适合精细化研究氧化还原前锋的动力学过程、矿物沉淀对渗透率的反馈效应
2.3 计算量的”三重爆炸”
反应-输运模拟的计算量来自三个维度的乘积：
空间维度：三维盆地尺度模型，网格数从几万到几百万不等。鄂尔多斯盆地铀矿勘查中，一个中等分辨率的区块模型（10km×10km×500m）采用100m×100m×10m的网格，网格数约50万。如果要刻画氧化还原前锋的精细结构（前锋宽度可能只有几米到几十米），需要局部加密到米级，网格数轻松破千万。
化学维度：铀矿成矿涉及的化学组分包括：主要离子（Na⁺、Ca²⁺、Mg²⁺、Cl⁻、SO₄²⁻、HCO₃⁻）、铀物种（UO₂²⁺、UO₂(CO₃)₂²⁻、UO₂(CO₃)₃⁴⁻等，取决于pH和碳酸根浓度）、还原剂（Fe²⁺、H₂S、有机质）、矿物相（沥青铀矿、黄铁矿、赤铁矿、方解石、石英、黏土矿物）、气体（O₂、CO₂、CH₄、H₂S）。一个中等复杂度的模型可能追踪20-50个组分、10-30个矿物相。
时间维度：砂岩型铀矿成矿时间尺度是百万年级别。模拟步长受CFL条件限制（dt < dx/v，v是孔隙流速）。在渗透性较好的砂岩中，地下水流速可能达到每年几米到几十米，100m网格的CFL步长约1-10年。百万年模拟需要10万-100万个时间步。如果采用隐式格式，步长可以放宽，但每个时间步的矩阵求解更耗时。
三个维度相乘：50万网格 × 30组分 × 10万步 = 1.5万亿次化学平衡计算。这是中等规模模型。如果是千万网格、50组分、隐式求解的大型三维模型，计算量再翻一个数量级。
三、热力学平衡计算：每一个时间步都在解一个非线性方程组
反应-输运模拟中最耗时的子任务，是化学平衡计算。在每个网格、每个时间步，都需要求解一组非线性方程，确定所有水溶液物种的活度、矿物的饱和状态、气体的逸度。
3.1 两种数学策略
目前有两种主流的热力学平衡算法：
质量作用定律法（LMA, Law of Mass Action）：以主要物种（primary species）为未知数，次生物种（secondary species）通过质量作用方程表达为主要物种的函数。求解质量平衡方程（每个元素的总量守恒）和电荷平衡方程。PHREEQC、MINTEQA2采用这种方法。优点是方程数量少（未知数=主要物种数），但质量作用方程的非线性较强，Newton-Raphson迭代的雅可比矩阵条件数可能很差。
吉布斯自由能最小化法（GEM, Gibbs Energy Minimization）：直接最小化系统的总吉布斯自由能，约束条件是元素总量守恒和非负摩尔数。GEMS、Perple_X、FactSage采用这种方法。优点是可以自然处理多相平衡（矿物、熔体、流体共存），避免了LMA法中”选择哪些矿物参与平衡”的主观性。缺点是优化问题的维度更高（变量数=所有物种数+所有矿物相数），对内存需求更大。
2024年玉龙铜矿水-岩反应热力学平衡模拟的论文，同时使用了LMA和GEM两种方法做对比验证。结论是：对于简单体系（少量矿物相），两种方法结果一致；对于复杂体系（多固溶体、多矿物相竞争），GEM更稳健，但计算时间增加30%-50%。
3.2 数值求解的痛点
热力学平衡计算的Newton-Raphson迭代，面临几个经典难题：
收敛性问题：当溶液接近矿物饱和边界（饱和指数SI≈0）时，矿物沉淀量对SI极其敏感，雅可比矩阵病态，迭代可能振荡或发散。PHREEQC通过”打靶法”（pivoting strategy）和”欠松弛”（under-relaxation）来缓解，但复杂模型仍然可能不收敛。
活度系数计算：高离子强度溶液（如卤水、尾矿库渗滤液）中，Debye-Hückel模型失效，需要Pitzer模型或B-dot模型。Pitzer模型涉及二元、三元相互作用参数，计算量比Debye-Hückel高一个数量级。
表面络合模型：矿物表面的质子化/去质子化、离子吸附，需要额外求解表面电荷平衡和电位（扩散双层模型）或多位点络合（CD-MUSIC模型）。这些非理想效应让Newton-Raphson的收敛域更窄。
在三维反应-输运模拟中，这些痛点被空间维度放大：如果一个网格的化学平衡不收敛，整个时间步可能回退重来；如果10%的网格频繁不收敛，总计算时间可能翻倍。
四、蒙特卡洛不确定性分析：成矿预测的概率量化
反应-输运模拟的结果，对输入参数极其敏感。渗透系数的不确定性可能来自实验室测定的误差（一个数量级）、地层非均质性（孔隙度从5%到30%）、还原剂含量（有机质丰度空间变异）。直接取参数中值做”最佳估计”模拟，得到的矿体位置可能偏离真实位置数百米。
这就是为什么蒙特卡洛不确定性分析成为成矿预测的标准流程。
4.1 参数敏感性分析
第一步是识别关键参数。通过拉丁超立方采样（LHS）或Sobol序列，在参数空间内生成数百到数千组样本，每组参数跑一次正演模拟。然后分析输出（如矿体厚度、铀品位、前锋线位置）对各个输入参数的响应。
砂岩型铀矿模拟中，最敏感的参数通常是：
- 渗透系数（控制地下水流速和前锋推进速度）
- 还原剂容量（控制氧化还原前锋的滞留时间）
- 铀源浓度（控制矿体品位）
- 孔隙度（控制溶质运移的储集空间）
- 弥散度（控制前锋的弥散宽度）
成都理工大学夏彧（2016）的博士论文《层间氧化带砂岩型铀矿成矿条件的模拟研究》中，作者用蒙特卡洛方法测试了渗透系数和孔隙度的敏感性，发现渗透系数的变化对前锋推进距离的影响是孔隙度的3-5倍。
4.2 概率预测与风险量化
第二步是用蒙特卡洛结果构建概率分布。例如，跑1000次模拟后，可以得到：
- 矿体存在概率图（每个网格出现矿化的频率）
- 铀品位累积分布函数（P10/P50/P90分位数）
- 前锋线位置的置信区间（95%置信度下的前锋线包络）
这对矿产勘查决策至关重要。如果一个区块的矿化概率<20%，即使”最佳估计”显示有矿，也不值得投入钻探。反之，如果概率>60%，即使最佳估计品位一般，也可能值得进一步验证。
4.3 蒙特卡洛的算力需求
蒙特卡洛的特点是embarrassingly parallel but computationally expensive。每次模拟完全独立，天然适合并行——1000次模拟可以在100核上同时跑，总时间压缩到串行的1/100。但每次模拟本身可能就需要几小时到几天（三维反应-输运），1000次的总量是惊人的。
以中等规模三维TOUGHREACT模拟为例：
- 单次模拟：50万网格 × 30组分 × 10万步，在32核服务器上跑约3-7天
- 1000次蒙特卡洛：串行需要3000-7000天（约8-19年）
- 并行到100核：压缩到30-70天（约1-2个月）
- 并行到500核：压缩到6-14天
这就是为什么大型地质调查机构（如美国USGS、中国核工业北京地质研究院）需要集群级算力来支撑成矿预测。
五、同位素混合模型：追踪地下水的”指纹”
反应-输运模拟告诉我们”发生了什么”，但不一定告诉我们”为什么发生”。要验证模型的可靠性，需要用同位素示踪来约束边界条件和反应路径。
5.1 同位素分馏的基本原理
同位素是同一元素的不同核素（质子数相同、中子数不同）。由于质量差异，同位素在物理化学过程中会发生分馏：
- 平衡分馏：热力学平衡状态下，轻重同位素在不同相之间的分配差异。例如，¹⁸O在水和方解石之间的平衡分馏系数与温度相关，是著名的氧同位素地质温度计。
- 动力学分馏：非平衡过程（如扩散、蒸发、生物代谢）中，轻同位素优先参与快速反应，导致分馏。例如，细菌还原硫酸盐时优先利用³²S，剩余硫酸盐富集³⁴S。
对于砂岩型铀矿，最有价值的同位素体系包括：
- 铀同位素（²³⁸U/²³⁵U）：²³⁸U优先被还原，还原沉淀的沥青铀矿相对富集²³⁸U，可用于识别还原环境
- 硫同位素（³⁴S/³²S）：硫酸盐还原过程中³⁴S分馏显著，可用于示踪还原剂的来源（生物成因vs热液成因）
- 碳同位素（¹³C/¹²C）：有机质氧化和碳酸盐溶解的¹³C特征不同，可用于区分碳来源
- 锶同位素（⁸⁷Sr/⁸⁶Sr）：示踪流体-岩石相互作用的时限和路径
5.2 同位素混合模型
在地下水系统中，不同端元（如大气降水、古湖水、深部卤水、油气流体）具有不同的同位素组成。混合模型通过质量平衡方程，计算混合比例：
δ_mix = f₁δ₁ + f₂δ₂ + … + fₙδₙ
其中fᵢ是第i个端元的体积分数，δᵢ是其同位素组成。对于二元混合，直接在δ-δ图解上画直线，样品落点指示混合比例。对于多元混合，需要结合其他地球化学指标（Cl⁻/Br⁻、Na⁺/Ca²⁺等）做联合约束。
2025年《中国环境科学》青年学术会议上，中国地质大学（武汉）刘鹏教授的报告”矿物吸附锑过程中的同位素分馏分子机制”，展示了第一性原理计算（DFT）如何预测同位素分馏系数——这是从宏观模拟向微观机理的跨越。
六、四级配置方案：从水文地球化学研究生到铀矿勘查国家队
结合反应-输运模拟、热力学平衡计算、蒙特卡洛不确定性分析的三重需求，以及UltraLAB产品目录，我梳理了四个层级的配置。
入门级：水文地球化学研究生（一维柱实验模拟、反向建模）
| 配件 | 推荐型号 | 说明 |
|---|---|---|
| CPU | Intel Ultra 285K（24核，P核最高5.5GHz） | PHREEQC批次反应主要串行，单核性能强不吃亏 |
| 内存 | 64GB DDR5-5600 | 支撑中等复杂度化学体系（20-30组分） |
| 存储 | 2TB NVMe SSD | PHREEQC输入输出文件快速读写 |
| GPU | RTX 5090 32G（标配） | 地球化学模拟纯CPU计算，GPU无加速；标配GPU可用于后期ML辅助参数优化 |
| 对应UltraLAB型号 | A350 15364-MBT（Intel Ultra 265K + 64GB + 2TB，标配双RTX 5090D 24G）或 A350 15696-MCE（Ultra 285K + 96GB + 4TB，标配RTX Pro5000 48G） | 入门级AI/图形工作站，标配GPU可作为通用算力储备 |
这套配置适合研究生学习PHREEQC基础操作、跑一维运移-反应模拟、做课程项目中的反向地球化学建模。能处理几百个网格的一维问题，但三维模型跑不了。
进阶级：课题组主力工作站（中等三维模型、参数敏感性扫描）
| 配件 | 推荐型号 | 说明 |
|---|---|---|
| CPU | AMD Ryzen Threadripper 7970X（32核64线程，基础4.0GHz/最高5.3GHz） | 多核并行参数扫描+蒙特卡洛 |
| 内存 | 128GB DDR5-5600 | 支撑中等三维模型（50万网格×30组分） |
| 存储 | 4TB NVMe SSD ×2 | TOUGHREACT输出文件动辄几十GB |
| GPU | 双RTX 4090 48G（标配） | 地球化学模拟无GPU加速；标配双GPU可用于后期ML辅助参数优化 |
| 对应UltraLAB型号 | AE450 145128-MC2T（AMD Threadripper 7970X + 双RTX 4090 48G水冷 + 128GB + 4TB） | 双GPU水冷图形工作站，适合课题组主力计算+后期ML扩展 |
这个级别可以支撑砂岩型铀矿区块尺度的三维反应-输运模拟（如10km×10km×500m的粗分辨率模型），以及数百组参数的敏感性分析。128GB内存是必要的，因为三维网格×组分×时间步的中间结果需要大量内存缓存。
专家级：成矿预测核心平台（高分辨率三维模型、全盆地蒙特卡洛）
| 配件 | 推荐型号 | 说明 |
|---|---|---|
| CPU | Xeon金牌6530 ×2（64核，2.7-4.0GHz） | 双路服务器级，MPI并行扩展到64核 |
| 内存 | 384GB DDR5-4800 | 高分辨率模型（百万网格）+大型化学体系 |
| 存储 | 8TB NVMe SSD + 20TB HDD阵列 | 海量模拟结果存档+地质数据库 |
| GPU | 不需要 | 除非后期用机器学习辅助参数优化，否则GPU非刚需 |
到了这个级别，你基本是在做鄂尔多斯盆地、伊犁盆地、松辽盆地级别的铀矿成矿预测。可以跑高分辨率三维模型（局部加密到米级）、全盆地流场模拟、数千次蒙特卡洛概率预测。384GB内存是必要的，因为TOUGHREACT的算子分裂求解需要存储巨大的中间状态数组和稀疏矩阵。
集群级：国家级铀矿勘查数据中心（超大规模蒙特卡洛、多盆地耦合）
| 配件 | 推荐型号 | 说明 |
|---|---|---|
| 计算节点 | GX668 227640-MD4E（2×Xeon金牌6530 + 640GB DDR5 + 4TB NVMe，标配4×RTX Pro6000D 84G） | 多节点MPI并行，每节点负责一组蒙特卡洛样本 |
| 存储系统 | 并行文件系统（Lustre/BeeGFS）+ PB级磁带库 | 全国铀矿勘查数据的长期存档 |
| 网络 | InfiniBand HDR 200Gbps | MPI节点间低延迟通信 |
| 总节点数 | 50-200节点 | 支撑数千次蒙特卡洛并行+多盆地同时模拟 |
| 对应UltraLAB型号 | GX668 2271T-MC8A（2×Xeon金牌6530 + 1TB DDR5 + 4TB NVMe，标配8×A100 80GB）或 GX668 227512-MC4A（4×A100 80GB） | 多GPU集群计算平台，标配GPU适合后期ML辅助参数优化+通用高并发计算 |
这就是中国核工业集团、美国能源部（DOE）级别的铀矿勘查计算平台了。美国怀俄明盆地、得克萨斯州南部的ISR（原位回收）铀矿，在环评阶段需要跑数千次蒙特卡洛来评估地下水污染风险——没有集群级算力，这活干不完。
七、避坑指南：反应-输运模拟研究者最容易踩的六个坑
坑一：热力学数据库不一致
地球化学模拟的结果，强烈依赖于热力学数据库（log K、活度系数参数、矿物焓/熵数据）。不同数据库（wateq4f.dat、llnl.dat、minteq.dat）对同一反应的log K可能相差0.5-2个单位，导致矿物饱和指数完全不同。
解决方案：文献调研你的温度-pH范围，选择专门校准过的数据库。高温地热系统用SPRONS或THERMOCALC数据库，低温地下水用wateq4f或llnl。跨数据库比较是审稿人的必考点，提前做。
坑二：算子分裂的耦合误差
PHREEQC和TOUGHREACT默认使用算子分裂（先运移后反应）。如果时间步长太大，运移和反应的非同步会引入显著误差——特别是在反应前锋附近，化学反应速率极高，分裂误差会让前锋位置偏移。
解决方案：用Courant数（Cr = v·dt/dx）检验，Cr < 0.5通常安全。对于强耦合问题（如矿物沉淀堵塞孔隙），改用CrunchFlow的全耦合求解，或TOUGHREACT的自适应时间步长。
坑三：忽略矿物动力学
很多初学者把所有矿物设为瞬时平衡。但现实中，矿物溶解/沉淀受动力学控制——石英的溶解速率极慢（10⁻¹⁴ mol/m²/s），方解石则快得多（10⁻⁶ mol/m²/s）。如果模型中让所有矿物瞬时平衡，会高估反应速率，预测的矿体位置可能提前。
解决方案：文献查阅关键矿物（沥青铀矿、黄铁矿、方解石）的动力学速率常数，用动力学表达式（Transition State Theory或empirical rate law）替代平衡假设。
坑四：蒙特卡洛采样不足
1000次蒙特卡洛听起来很多，但如果参数维度高（>10个敏感参数），1000次可能远远不够覆盖参数空间。Latin Hypercube采样在高维时仍然稀疏，Sobol序列更好但实现复杂。
解决方案：先用少量样本（100-200次）做初步敏感性分析，识别关键参数（通常3-5个），再针对关键参数做密集采样（1000-5000次）。高维问题考虑用代理模型（Kriging、高斯过程）替代正演模拟。
坑五：边界条件设置不合理
反应-输运模拟的结果对边界条件极其敏感。一个常见的错误是：把氧化还原前锋设置为固定边界（Dirichlet条件），而实际上前锋是模型内部自发演化的结果。
解决方案：氧化还原前锋应该由模型内部的水-岩反应自然产生，而不是人为设定。补给边界只设定水化学组成（如大气降水成分），让模型自己算出前锋位置和形态。
坑六：存储管理混乱
三维反应-输运模拟的输出文件体积极大。一个中等规模TOUGHREACT模拟，如果每个时间步都保存全部网格的化学状态，输出轻松超过100GB。如果没有自动化存储管理，硬盘很快满。
解决方案：只保存关键时间步（如每1000年或每10000步）的完整状态，中间步只保存监测点（如钻孔位置）的数据。用HDF5格式替代ASCII，体积缩小5-10倍。定期归档旧结果到HDD或磁带。
八、从同位素讲习班到数字孪生：中国地球化学的四十年算力长征
1983年北大那间教室里，郑淑慧副教授用粉笔在黑板上推导氢氧同位素分馏公式。学生们用手摇计算器算活度系数，用绘图纸画δ¹⁸O-δD图解。一个批次反应的平衡计算，可能需要一整天。
四十年后，同样的计算在一台笔记本上只需要毫秒级。但地球化学家的问题规模也指数级增长：从一维柱实验到三维盆地尺度，从单一组分到多相多组分耦合，从确定性预测到蒙特卡洛概率量化。
砂岩型铀矿的成矿预测，是中国能源安全的战略性需求。中国已探明的砂岩型铀矿主要分布于鄂尔多斯、伊犁、二连、松辽等盆地，但仍有大量空白区等待勘查。反应-输运模拟+蒙特卡洛不确定性分析+同位素示踪约束，这套方法论正在把”找矿靠运气”变成”找矿靠算力”。
1983年北大的讲习班，为中国地球化学研究奠定了方法论基础。今天的工作站和集群，则在用算力解码氧化还原前锋的推进轨迹，用概率量化矿体的藏身之处。从”铀随水去”到”水去铀留”，这场地下水里的氧化还原接力贯穿了地球化学演化的漫长历史。用计算机重演它，可能需要跑几天到几个月。但比起找到一个新的铀矿床，这点等待不算什么。
参考文献与延伸阅读：
- 聂逢君, 严兆彬, 夏菲, 何剑锋, 张成勇, 封志兵, 张鑫, 杨东光, 陈梦雅, 谈顺佳, 张进, 康世虎, 宁君, 杨建新, 申科峰, 蔡建芳（2025）. 砂岩型铀矿的”双阶段双模式”成矿作用. 沉积与特提斯地质, 45(2): 1-15.
- 金若时, 苗培森, 李建国, 林双幸, 王高尚, 蒋乐, 惠小朝, 余安（2023）. 中国北方砂岩型铀矿大规模成矿作用与成矿模型. 铀矿地质, 39(4): 193-210.
- 夏彧（2016）. 层间氧化带砂岩型铀矿成矿条件的模拟研究. 成都理工大学博士学位论文.
- 吴德海, 张金带, 李建国, 等（2022）. 砂岩型铀矿CO₂+O₂地浸采铀的反应运移数值模拟. 中国科学:技术科学, 52(8): 1456-1470.
- 刘武生, 牛霆, 贾立城, 杨喆, 李伟涛, 李真真, 史清平（2021）. 酸法原地浸出开采铀矿的水文地球化学模拟. 核工业北京地质研究院技术报告.
- Parkhurst D.L., Appelo C.A.J.（2013）. Description of input and examples for PHREEQC version 3. USGS Techniques and Methods, book 6, chap. A43.
- Xu T., Spycher N., Sonnenthal E., Zhang G., Zheng L., Pruess K.（2011）. TOUGHREACT User’s Guide: A Simulation Program for Non-isothermal Multiphase Reactive Geochemical Transport in Variably Saturated Geologic Media. LBNL-62720E.
- Steefel C.I., Appelo C.A.J., Arora B., et al.（2015）. Reactive transport codes for subsurface environmental simulation. Computational Geosciences, 19(3), 445-478.
- 张金带, 修高中, 陈安平, 王成（2004）. 斑岩铜矿蚀变分带与成矿机制：玉龙矿床水-岩反应热力学平衡模拟例析. 地质科学, 39(2): 196-208.
- 陈留勤, 王京平, 李遥, 何立羹, 黄天明, 皮坤福, 修伟, 张俊文, 苏正旗, Amjad Danyal, 胡晓农, 李俊霞, 王健, 王雨桐, 臧春艳, 钟承昊, 罗豪豪（2025）. 全流程PHREEQC水文地球化学模拟及高阶耦合技术实战. 培训课程讲义.
- 郭华明, 高金尉, 苏明, 周伟青, 陈林捷, 刘星, 朱振利（2024）. 地下水环境地球化学过程与污染修复. 科创中国 成果库.
- 中国矿物岩石地球化学学会（1983）. 稳定同位素地球化学基础理论讲习班. 中国矿物岩石地球化学学会年会资料.
