本文档报告了一项关于三维粗糙壁岩石裂隙中溶质运移预测的原创性研究。该研究由Na Huang, Zhenhui Wang, Yujing Jiang, Ju Wang, Bo Li等人完成，作者单位包括山东科技大学土木工程与建筑学院、长崎大学研究生院工程学研究科、核工业北京地质研究院以及同济大学土木工程学院。该研究发表于Journal of Hydrology期刊，于2026年5月23日在线发布。

学术背景 本研究属于水文地质学、岩石力学及多孔介质传输物理的交叉领域，核心关注裂隙岩体中的流体流动与溶质运移问题。在地壳中，裂隙是流体流动和溶质迁移的关键通道，尤其对于基质孔隙度和渗透率较低的结晶岩体而言。准确刻画岩石裂隙中的溶质运移过程，对于解决地下水污染修复、地热能开采和放射性废物地质处置等众多地学工程问题至关重要。

然而，由于单一裂隙固有的几何复杂性（如壁面粗糙度、接触点分布）及其导致的流动异质性，直接预测其传输过程极具挑战。传统的平滑平行板模型过于简化，无法捕捉实际粗糙裂隙中的非费克（non-Fickian）输运行为（如早期突破和长拖尾现象）。虽然已有多种复杂模型（如分数阶对流-弥散方程、连续时间随机游走、活动-非活动域模型等）被提出以描述此类异常输运，但这些模型通常需要拟合多个参数，且参数物理意义不明确，限制了其在实际应用中的预测能力。因此，开发一种能够基于可直接测量的几何参数、无需复杂反演即可直接预测溶质运移的模型，具有重要的科学意义和实际价值。

本研究旨在建立一个直接预测框架，通过明确纳入表征几何非均质性的两个关键参数——局部力学开度的相对标准偏差（Relative Standard Deviation, RSD）和接触面积比（Contact Area Ratio, CR），来预测层流条件下三维粗糙壁岩石裂隙中的溶质运移。其核心目标是建立一个半经验解析模型，该模型能够仅基于几何参数（RSD, CR）和水力开度，直接预测突破曲线（Breakthrough Curves, BTCs），从而将传统对流-弥散方程（ADE）从描述性工具转变为预测性工具。

详细研究流程 本研究包含一个系统且逻辑严密的工作流程，主要可分为以下几个步骤：

1. 三维粗糙壁裂隙模型的构建：

   • 研究对象与方法：研究构建了一系列具有不同几何特征的三维粗糙壁裂隙数值模型。裂隙表面采用分数布朗运动（Fractional Brownian Motion, FBM）方法生成，以捕捉天然裂隙形貌的统计自仿射特性和长程相关性。具体使用了Liu等人（2004）提出的增强型连续随机添加算法来合成具有不同赫斯特指数（H = 0.5–0.7）的裂隙表面，其对应的节理粗糙度系数（JRC）范围为3.18至15.17，覆盖了天然岩石裂隙的典型粗糙度。
   • 模型生成与参数化：通过控制上下裂隙面的侧向偏移和垂直分离，生成具有可变开度场的非匹配裂隙。通过改变垂直分离距离，使力学开度算术平均值（b）在0.2至1.4 mm范围内变化，同时将接触面积比（CR）控制在0.25以下，以避免流道堵塞。最终，从生成的裂隙中截取0.1 m × 0.05 m的截面作为流动与运移分析模型，共构建了从F1到F5的5组不同粗糙度的模型系列，每个系列包含多个不同平均开度和CR的样本。研究提取了每个裂隙模型的两个关键几何参数：RSD（仅基于空隙空间计算）和CR。所有模型的几何特性详见表1。

2. 流动与溶质运移的数值模拟：

   • 控制方程与数值方法：采用商业有限元软件COMSOL Multiphysics求解控制方程。稳态单相流动由Navier-Stokes方程和连续性方程描述。溶质运移通过求解对流-扩散方程进行模拟。流体性质采用水（密度1000 kg/m³，动力粘度0.001 Pa·s），溶质采用典型保守溶质（如Cl⁻，分子扩散系数D0 = 2.03 × 10⁻⁹ m²/s）。
   • 边界条件与参数设置：入口设定恒定体积流量（Q），出口压力设为零以驱动流动。裂隙壁面及接触区域设置为无滑移边界。通过控制Q，将雷诺数（Re）保持在0.01，确保流动处于线性层流状态，排除惯性效应影响。Péclet数（Pe）为4.93，低于渐近输运 regime 的临界值，因此观测到的非费克输运行为可归因于开度非均质性，而非裂隙长度不足。溶质注入在流场稳定后开始，入口为恒定浓度边界（C0 = 1），初始浓度为零。
   • 网格划分与验证：裂隙域被离散为约2.0×10⁷个非结构四面体网格。在裂隙壁面和接触区域附近采用精细网格（0.03 mm），内部区域采用较粗网格（0.15 mm）。通过网格敏感性分析确认了网格分辨率足够。设定0.01 mm的阈值开度，低于此值的局部开度被视为接触区域以简化网格。该数值模拟方法通过透明树脂3D打印制作的物理模型实验进行了验证，对比显示数值模拟在溶质羽流形态和突破曲线关键特征（突破时间、峰值浓度）上与实验结果高度吻合，证实了数值模型的有效性（图4）。

3. 突破曲线分析与经典模型拟合：

   • 数据处理：基于数值模拟得到的出口断面流速和浓度分布，通过积分计算出口溶质通量，进而得到通量加权浓度，并归一化为C‘和孔隙体积（PV），从而获得所有裂隙模型的BTCs（图6）。
   • 反演分析：使用经典的一维对流-弥散方程（ADE）及其解析解（Ogata and Banks, 1961）对每条BTC进行反演拟合，以获取两个关键运移系数：运移速度（v）和水动力弥散系数（D）。拟合优度通过均方根误差（RMSE）评估。所有模型的RMSE均低于0.03，表明ADE模型能够高保真地再现所有BTCs，有效捕捉了其非费克输运特征。此步骤是模型校准，为后续的参数化建立基础。

4. 运移系数与几何参数的参数化关系建立：

   • 数据归一化：将反演得到的v和D分别除以由平滑平行板模型理论公式计算出的对应值（vp和dp），得到归一化运移速度（v/vp）和归一化弥散系数（d/dp）。vp和dp根据立方定律和泰勒弥散理论公式计算（公式9，10）。
   • 单参数关系分析：首先分别分析v/vp和d/dp与RSD、CR的单变量关系（图7，图8）。结果表明：
     • v/vp随RSD和CR的增加而单调递减，反映了几何非均质性（粗糙度和接触）导致的流动曲折度和通道化加剧，降低了宏观运移速度。与RSD呈线性关系（v/vp = 1 – 0.5855×RSD），与CR呈幂函数关系（v/vp = 1 – 0.8770×CR^0.3479）。
     • d/dp随RSD和CR的增加而单调递增，范围在1到140之间，表明几何非均质性显著增强了机械弥散。与RSD呈指数关系（d/dp = e^(5.8123×RSD)），与CR呈幂函数关系（d/dp = (1 + 60.0578×CR)^1.9258）。
   • 双参数综合模型建立：基于单参数分析，建立了同时包含RSD和CR的综合性参数化函数（公式20，21）：
     • f1(RSD, CR) = v/vp = 1 – λ×RSD – β×CR^θ （λ = 0.4536, β = 0.3272, θ = 0.4827, R² = 0.9789）
     • f2(RSD, CR) = d/dp = e^(χ×RSD) × (1 + κ×CR)^ω （χ = 4.3110, κ = 100, ω = 0.4131, R² = 0.8432） 这两个函数满足物理约束：当RSD→0且CR→0（平滑平行板）时，函数值趋近于1；且随RSD和CR增加，f1单调递减，f2单调递增。

5. 预测模型的集成与验证：

   • 模型集成：将上述参数化函数f1和f2代入经典的ADE解析解（公式8），得到一个修正后的解析解（公式14）。该模型仅需要输入几何参数（RSD, CR）和水力开度（bh，可通过已有经验公式由力学开度、JRC等估算）以及水力梯度，即可直接预测BTC，无需进行详细的流场和浓度场计算或BTC反演拟合。
   • 内部验证：使用该模型对研究中所有裂隙模型的BTC进行预测，并与反演拟合结果比较。预测的RMSE值略高于反演拟合，但仍全部低于0.04，表明模型具有可靠的预测精度（图10）。
   • 外部验证与对比：为了评估模型在更广泛裂隙上的适用性，研究选取了Li等人（2008）实验中的三个剪切裂隙样本（J1, J2, J3）进行数值模拟生成基准BTC，并将本模型预测结果与现有其他模型（Zhou et al., 2023; Dou et al., 2018; Wang et al., 2012）的预测进行对比（图12）。结果显示，本模型预测的BTC与基准数据吻合良好（RMSE < 0.043），而其他模型由于未考虑接触面积（Zhou et al., 2023）、基于二维剖面（Dou et al., 2018）或主要忽略了几何非均质性（Wang et al., 2012），预测结果出现不同程度的偏差。这凸显了本模型在同时考虑开度变异性和接触面积方面的优势。

主要研究结果 1. 几何非均质性是控制溶质运移的关键因素：在层流条件下，观测到的非费克输运行为（早期突破和长拖尾）主要归因于裂隙的几何非均质性，而非惯性涡流。粗糙度和接触面积共同塑造了高度异质性的流场（形成优势流道和滞流区）和浓度场，从而导致异常运移。 2. 运移参数与几何参数存在明确的定量关系：研究成功地将归一化运移速度（v/vp）和归一化弥散系数（d/dp）与两个可测量的无量纲几何参数（RSD和CR）进行了参数化。v/vp随RSD和CR增加而减小，d/dp随RSD和CR增加而增大。这为基于几何特征预测运移行为提供了直接、定量的桥梁。 3. 建立了半经验解析预测模型：通过集成上述参数化关系，本研究提出了一个修正的ADE框架（公式14）。该模型是一个半经验解析模型，其核心创新在于将几何异质性通过f1和f2两个校正函数直接嵌入经典解析解中。 4. 模型具有高精度和广泛适用性：模型在内部数据集和外部独立实验样本上均表现出良好的预测性能。它能够准确再现不同粗糙度、开度和接触面积比的裂隙中的溶质突破曲线。更重要的是，当几何非均质性消失（RSD→0, CR→0）时，模型能自适应地退化为经典的平滑平行板模型解，体现了其理论的完备性。 5. 接触面积的作用得到量化：与以往忽略接触或仅基于二维/无接触模型的研究相比，本研究明确量化了接触面积比（CR）对运移的显著影响。接触区域不仅进一步降低了有效运移速度，而且极大地增强了溶质的弥散（d/dp可高达平滑模型的140倍），这是先前模型未能充分捕捉的。

结论与意义 本研究成功开发了一个直接预测三维粗糙壁岩石裂隙中溶质运移的框架。该框架通过明确引入并参数化相对标准偏差（RSD）和接触面积比（CR）这两个几何非均质性指标，建立了一个半经验解析模型。该模型能够仅基于可测量的几何参数和水力条件，直接、准确地预测层流条件下的溶质突破曲线，无需进行复杂的数值模拟或反演拟合。

其科学价值在于：1）深化了对裂隙几何结构如何控制溶质运移机理的理解，特别是量化了开度变异性和接触面积的独立与耦合效应；2）提供了一种将复杂三维几何信息与宏观运移参数联系起来的参数化方法，为从结构预测功能提供了新途径；3）所提模型兼具物理明确性和预测实用性，是对现有复杂非费克输运模型集的一个重要补充。

其应用价值在于：该模型为地下水污染风险评估、地热资源开发、核废料处置库性能评估等领域的工程实践提供了一个高效、实用的预测工具。现场可通过钻孔成像、激光扫描等技术获取裂隙表面形貌和接触信息，进而估算RSD和CR，从而直接预测溶质运移行为，降低了对大量水力试验和复杂数值模拟的依赖。

研究亮点 1. 研究对象的全面性：首次在预测模型中同时系统地纳入了裂隙开度变异性（通过RSD表征）和接触面积（通过CR表征）这两个关键几何控制因素，更贴近实际裂隙（尤其是受应力作用发生闭合或剪切位移的裂隙）的几何特征。 2. 方法论的创新性：创造性地将几何参数与经典ADE模型的运移系数通过归一化和参数化函数直接关联，构建了一个“基于几何的预测”框架。该方法避免了多参数反演的非唯一性问题，参数物理意义清晰。 3. 模型的预测性与普适性：所开发的模型是解析/半解析形式的，具有直接预测能力，且通过严格的数值实验和外部数据验证，证明了其在较宽几何参数范围内的有效性和优于部分现有模型的性能。 4. 清晰的物理图像与退化一致性：模型具有明确的物理基础（校正函数f1, f2的单调性符合物理直觉），并且能平滑退化为经典理论解，体现了模型在理论上的严谨性。

其他有价值内容 研究还讨论了模型的局限性，包括：1）目前仅适用于层流（线性达西流）条件，在更高雷诺数（非达西流或湍流）下的性能有待验证；2）CR是一个全局参数，未能反映接触区域的空间分布特征，未来需探索能描述接触空间结构的稳定参数；3）模型针对单向流、连续注入和单一孤立裂隙校准，其在径向流、脉冲注入或裂隙网络中的直接适用性尚未验证；4）模型假设溶质为保守性，未考虑吸附、溶解、沉淀等反应性输运过程。这些讨论为未来研究指明了方向。