(a)裂缝及周围损伤区; (b)损伤区岩石中裂缝分布情况; (c)岩石中裂缝密度的变化规律
图1 裂缝及断层损伤区的结构示意图(修改自文献[13])
Fig. 1 Structure of fracture and fault damage zone (modified from Ref. [13])
北京大学学报(自然科学版) 第62卷 第1期 2026年1月
Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 62, No. 1 (Jan. 2026)
doi: 10.13209/j.0479-8023.2025.027
国家自然科学基金(42372305)和国家重点研发计划(2021YFA0716000)资助
收稿日期: 2025–01–23;
修回日期: 2025–04–16
摘要 根据现场示踪剂试验的统计结果, 裂隙介质中的有效基质扩散系数会随着观测尺度的增加而增加。为了研究其背后的机制, 建立考虑裂缝周围损伤区的三维裂缝–基质模型, 利用溶质运移过程的数值模拟结果, 分析损伤区厚度、渗透率以及孔隙度对有效基质扩散系数的影响。结果表明, 当损伤区具有较大的厚度和孔隙度时, 会使得更多溶质进入岩石基质, 导致测得的有效基质扩散系数增大。对于天然裂缝来说, 损伤区内部岩石的渗透率及孔隙度会随着与裂缝面距离的减小而增大。当裂缝尺度增大时, 损伤区的厚度增加, 岩石的渗透率和孔隙度也整体上增大。不同观测尺度下的模拟结果表明, 有效基质扩散系数的尺度因子与观测尺度间存在正向的关系。该研究结果可以解释有效基质扩散系数尺度效应的部分机制, 有助于了解复杂裂隙–基质系统中的溶质运移过程。
关键词 基质扩散; 尺度效应; 损伤区; 溶质运移; 裂隙介质
在裂隙介质的溶质运移过程中, 基质扩散是重要的物理机制, 关于基质扩散的研究在污染物运移、核废料存储及能源勘探开发等领域具有重要意义[1]。有效基质扩散系数是反映岩石中基质扩散能力的关键参数, 已在大量现场试验及实验室研究中被测量[2–4]。然而, 现场示踪剂试验测得的有效基质扩散系数显著大于实验室测量结果, 在一些试验中, 二者的差异甚至达到几个数量级[5–6]。
近年来, 大量研究致力于解释产生有效基质扩散系数尺度效应的物理机制[7–10], 但其结果均无法与现场试验数据完全匹配。对此, 学术界存在两种主要观点。一部分学者认为有效基质扩散系数的尺度效应与增强的扩散作用有关, 如 Neretnieks[7]认为裂缝中存在的死水区延缓溶质的流动, 增强了溶质的分子扩散过程, Dai 等[8]认为岩石基质的非均质性使得计算得出的有效基质扩散系数大于平均值。另一部分学者认为增强的基质扩散现象在本质上与溶质到达的延缓作用相同, 如 Lei 等[9]提出粗糙的裂缝表面增加了溶质运移的路径, Liu 等[10]认为裂缝周围密集存在的小尺度裂缝增加了溶质运移路径的长度, 这些情况都会增加溶质运移的时间, 使得计算得出的有效基质扩散系数增大。
实际上, 对地质体中的天然裂缝及断层来说, 其周围一般存在一个有大量次级裂缝的区域[11–12], 即裂缝损伤区。在复杂的地质作用下, 该区域岩石发生破裂和变形, 形成大量密集且互相连通的裂缝, 其结构如图 1 所示, 可见损伤区岩石的裂缝密度明显大于远离裂缝面处未受破环的岩石。研究表明, 损伤区的裂缝密度会根据裂缝的尺度差异而发生改变, 且这些密集的小尺度裂缝改变了岩石的物理性质[11,13–14], 可能对溶质运移过程产生影响。Liu 等[10]研究了裂缝周围损伤区对溶质运移的影响, 但他们只考虑小尺度裂缝分形特性对溶质运移路径的影响, 忽略了损伤区内密集裂缝对岩石性质的影响, 也没有考虑不同尺度下裂缝损伤区内岩石性质的差异。
(a)裂缝及周围损伤区; (b)损伤区岩石中裂缝分布情况; (c)岩石中裂缝密度的变化规律
图1 裂缝及断层损伤区的结构示意图(修改自文献[13])
Fig. 1 Structure of fracture and fault damage zone (modified from Ref. [13])
本研究通过建立具有裂缝周围损伤区的三维裂缝–基质模型, 对其中的溶质运移过程进行数值模拟。为了研究损伤区对有效基质扩散系数的影响, 设置裂缝周围具有不同厚度及孔隙度的均匀损伤区模型, 用于分析当损伤区孔隙度和厚度改变时有效基质扩散系数的变化情况。为了模拟真实裂缝损伤区中溶质运移过程, 设置损伤区内部岩石的渗透率及孔隙度随与裂缝面距离变化的不均匀损伤区模型。最后, 模拟不同裂缝尺度下的溶质运移过程, 当裂缝尺度增加时, 设置的损伤区厚度、渗透率和孔隙度均出现增大趋势, 通过计算不同尺度下有效基质扩散系数的尺度因子, 分析产生有效基质扩散系数尺度效应的可能机制。
根据裂缝及断层损伤区的结构, 我们设置具有损伤区的三维裂缝–基质模型(图 2), 并对其中的溶质运移过程进行数值模拟。该模型为一个边长为400m 的立方体, 模型中间位置存在一条水平的光滑裂缝, 在其表面设置控制流体注入和抽出的注入井及监测井。由于天然裂缝周围存在具有较高密度
图2 考虑损伤区的裂缝–基质系统三维模型
Fig. 2 3D model of the fracture-matrix system with damage zone
裂缝的损伤区, 故在裂缝两侧设置对称分布且与裂缝面平行的损伤区岩石, 单侧岩石的厚度为 W。在损伤区之外为未受破坏的岩石, 其渗透率与孔隙度不随空间位置变化。
在实际的示踪剂试验中, 通常使用监测到的示踪剂浓度曲线来反演有效基质扩散系数。本研究根据数值模拟结果, 采用抽–注水井条件下溶质运移模型的解析方法[15]对溶质浓度曲线进行反演, 得到有效基质扩散系数。
1.2.1 抽–注水模型的解析解
根据 Webster 等[15]的理论研究, 在抽–注水井的溶质运移模型中, 流体的速度与空间位置有关, 形成一种双偶极子形状的流场, 在不考虑流线间溶质的相互交换及扩散时, 可以将溶质的流动过程近似地等效为多条单一流线的数值叠加, 监测点处的浓度曲线可表示为
, (1)
其中, n表示计算过程中叠加的流线数目(为了减少数值计算导致的误差, 设置 n=50), t表示时间, 
表示单一流线在监测点处的溶质浓度曲线, 其计算公式修改自 Tang 等[16]提出的简单二维解析模型表达式, 即
(2)
其中,
(3)
t0表示溶质注入的时间, C0表示溶质的浓度, u表示流体在裂缝中的速度, x表示单一流线的长度, Df为裂缝中的弥散系数, θ表示基质的孔隙度, Dm表示基质扩散系数, b表示裂缝开度的一半。
1.2.2 反演方法
在计算有效基质扩散系数时, 采用的反演算法为 Levenberg-Marquardt 方法。Levenberg-Marquardt是一种非线性的最小二乘方法[17], 可以在多次迭代并更新有效基质扩散系数的值后, 使得模拟结果与解析解得到的浓度曲线间误差的平方和最小。该方法通过引入一个调整因子, 在每一次的迭代过程中确定下一步的搜索方向并调整步长, 实现步长的自适应性。Levenberg-Marquardt 方法结合高斯–牛顿法和梯度下降法的优点, 能够快速收敛到最优解, 也能提高算法的稳定性。在本文的反演计算中, 除有效基质扩散系数外, 其他参数均为已知参数, 其中岩石孔隙度的设置与远离裂缝表面的未破坏岩石相同。
1.2.3 尺度因子
计算得出的有效基质扩散系数是一个绝对参数, 其值会受到示踪剂种类的影响[5]。为了明确造成有效基质扩散系数尺度效应的物理机制, 本文采用与 Zhou 等[5]类似的处理方法, 将尺度因子 F(反演得到的有效基质扩散系数与设置的基质扩散系数之比)作为表征基质扩散影响的参数, 其表达式为
, (4)
其中, 
为反演得出的有效基质扩散系数; 
表示设置的有效基质扩散系数, 也可以表示为实验室测得的基质扩散系数。
1.3.1 岩石性质均匀的损伤区模型
为了确定损伤区厚度和孔隙度对有效基质扩散系数的影响, 我们建立图 2 所示三维裂缝–基质模型。设置损伤区内部岩石性质不随空间位置变化, 并根据大量野外岩石样本的统计结果[5,14,18], 选取具有不同损伤区厚度 W (0.0001, 0.001, 0.01, 0.1 和 1m)和不同孔隙度 θ(0.001, 0.005, 0.01, 0.05 和 0.1)的情况进行数值模拟。模型参数设置见表 1。
表1 裂缝–基质模型的参数设置
Table 1 Parameters of the fracture-matrix model
参数名称参数值 井间距x (m)10 注入速率Q ( m3/s)0.0005 注入时间T0 (h)5 未破坏岩石孔隙度θ0.001 裂缝开度2b (m)0.002 裂缝渗透率kf (m2)8.33×10−18 基质扩散系数Dm (m2/s)1×10−8 基质渗透率km (m2)1×10−18
1.3.2 岩石性质随与裂缝面距离变化的不均匀损伤区模型
实际裂缝周围损伤区岩石的厚度、渗透率及孔隙度往往与裂缝的尺度有关。为了精确地描述不同尺度裂缝条件下损伤区岩石性质的变化情况, 并分析其对有效基质扩散系数的影响, 我们参考有关裂缝周围损伤区岩石性质变化规律的研究成果[14,19–20]。根据 Mitchell 等[14]对不同研究区地质体中损伤区裂缝分布情况的调查和统计结果, 损伤区岩石的渗透率变化规律可以表示为
, (5)
其中, k 表示损伤区岩石的渗透率; c 表示裂缝的密度; g=1.4832, 为根据统计数据得到的拟合参数; F0为损伤区在裂缝边界处的最大裂缝密度, 与岩石的性质有关; b=147.16, 为根据统计数据得到的拟合参数; D表示裂缝的位移量, 与裂缝规模存在正向的线性关系[19–20]; 
表示与裂缝面的距离。
根据上述不同尺度下损伤区岩石渗透率的研究结果, 我们设置图 3(a)所示不同尺度裂缝的损伤区渗透率模型, 其中 L表示裂缝的尺度, 远离裂缝面的未破坏岩石的渗透率为 1×10−18m2。可以发现, 随着与裂缝面距离减小, 损伤区渗透率逐渐增大, 在靠近裂缝面处达到最大值。
裂隙介质中岩石渗透率 k 与孔隙度 θ 之间存在如下对应关系[18]:
。 (6)
我们利用式(6), 设置图 3(b)所示不同尺度裂缝的损伤区孔隙度模型, 未破坏岩石的孔隙度为 0.0033。与损伤区渗透率的变化情况相同, 损伤区岩石孔隙度随着与裂缝表面距离的减小而增大。此外, 随着裂缝尺度增加, 损伤区的厚度、渗透率及孔隙度均呈上升趋势。
下面分析损伤区内部岩石为均匀介质时, 损伤区厚度及孔隙度对溶质运移过程的影响。在监测井处得到的随时间变化的溶质浓度曲线如图 4 所示, 可见随着损伤区孔隙度和厚度增大, 溶质浓度曲线的峰值逐渐降低, 表明岩石基质中发生更显著的分子扩散。图 4(b)显示, 当损伤区厚度为 0.0001 和0.001m 时, 曲线间存在非常小的差异, 此时基质扩散对损伤区厚度的变化不敏感; 当损伤区厚度为0.1 和 1m 时, 曲线间发生重合, 说明此时损伤区厚度对基质扩散的影响到达上限。
通过反演不同损伤区厚度及孔隙度情况下的溶质浓度曲线, 我们计算得出有效基质扩散系数的尺度因子随裂缝周围损伤区厚度及孔隙度增长的变化情况(图 5)。可以发现, 当损伤区孔隙度为 0.001 时, 由于损伤区与远离裂缝面的原始岩石具有相同的渗透率和孔隙度, 计算得出的尺度因子为 1, 即损伤区不对有效基质扩散系数产生影响。当损伤区厚度固定时, 随着孔隙度增大, 计算得出的尺度因子也呈现上升趋势。同时, 对于固定的损伤区孔隙度, 计算得出的尺度因子也随着损伤区厚度的增大而增大。这是因为增大的损伤区厚度及孔隙度为分子扩散提供了更大的空间, 允许更多的溶质进入岩石基质并发生停留, 从而使得岩石中发生的分子扩散过程更加显著, 计算得出的有效基质扩散系数也更大。此外, 当损伤区厚度为 0.0001m 时, 尺度因子随着损伤区孔隙度的变化并不显著, 说明此时岩石基质中分子扩散的强度对岩石孔隙度的变化不敏感。随着损伤区厚度增加, 尺度因子随损伤区孔隙度的变化趋势逐渐趋近于线性增长。当损伤区厚度为 1 和 0.1m 时, 两条曲线重合, 说明在 0.01～0.1m之间存在一个影响基质扩散的损伤区厚度上限。此时, 损伤区厚度的增加无法使更多的溶质进入基质并发生分子扩散, 这与观察溶质浓度曲线(图 4)得到的结论一致。
天然裂缝周围损伤区的渗透率及孔隙度会随着与裂缝面距离的减小而增大。我们模拟这种不均匀损伤区中溶质的运移过程, 设置的裂缝尺度为 10m, 损伤区厚度为 0.2m, 远离裂缝面未破坏岩石的渗透率为 1×10−18m2, 孔隙度为 0.0033, 靠近裂缝面岩石的渗透率为 1×10−16m2, 孔隙度为 0.0147。
图 6 显示, 当损伤区内岩石渗透率及孔隙度的变化符合天然裂缝周围损伤区的变化规律时, 其溶质浓度曲线与不考虑损伤区的结果有明显的差异。考虑损伤区岩石性质变化的溶质浓度曲线具有更低的峰值, 说明此时有更多的溶质进入岩石中, 并发生更显著的分子扩散。同时, 与无损伤区的溶质浓度曲线相比, 考虑损伤区时的溶质浓度曲线峰值降低 18%。对该曲线进行反演, 得出的有效基质扩散系数为 4.4×10−8m2, 尺度因子为 4.4, 有效基质扩散系数明显增大。上述结果说明, 这种岩石性质不均匀的损伤区会对溶质运移过程产生影响, 增大有效基质扩散系数的尺度因子。
图3 不同裂缝尺度下损伤区渗透率(a)及孔隙度(b)的变化
Fig. 3 Permeability (a) and porosity (b) of the damage zone for the different fracture scale
图4 损伤区孔隙度(a)及厚度(b)对溶质浓度的影响
Fig. 4 Effect of the damage zone porosity (a) and width (b) on the solute concentration
图5 损伤区厚度及孔隙度对尺度因子的影响
Fig. 5 Effect of the damage zone width and porosity on the scale factor
图6 不均匀损伤区模型的溶质浓度曲线
Fig. 6 Solute transport curve of the model with inhomogeneous damage zone
现场示踪剂试验结果表明, 有效基质扩散系数随着观测尺度的增大而增大。为了与不同观测尺度下的示踪剂试验结果进行对比, 我们设置不同的井间距来表示观测尺度的差异。通常情况下, 地下裂缝的尺度都等于或大于开展试验的井间距。为了精确地研究裂缝尺度对有效基质扩散系数的影响, 本文设置裂缝尺度与井间距相同。由于损伤区的厚度、岩石渗透率和孔隙度随裂缝尺度的增加而增加, 我们设置不同裂缝尺度下具有不同损伤区厚度的模型, 其损伤区内岩石渗透率及孔隙度随着与裂缝面的距离而变化。
对裂缝–基质系统中的溶质运移过程进行模拟, 计算得出的尺度因子随观测尺度变化的结果如图 7所示。可以发现, 当观测尺度小于等于 1m 时, 现场试验和模拟结果的尺度因子都趋近 1。这是因为此时损伤区的厚度小于 0.02m, 且损伤区岩石的孔隙度也与远离裂缝面的岩石间存在非常小的差异, 损伤区对有效基质扩散系数几乎没有影响, 模拟结果与不考虑损伤区的实验室试验结果相似。随着观测尺度增大, 模拟结果与现场试验结果的尺度因子均出现增长趋势。对于损伤区模型, 随着裂缝尺度增加, 损伤区厚度和孔隙度均呈现增加的趋势, 使得进入岩石的溶质增多, 在岩石中发生更显著的分子扩散, 导致计算得出的有效基质扩散系数增大。
与现场试验数据的对比结果显示, 在双对数坐标下的线性拟合结果中, 模拟数据具有 0.29 的斜率, 与现场试验数据 0.57 的斜率存在差异。一方面, 由于损伤区厚度对有效基质扩散系数的影响存在上限, 当观测尺度增大到一定程度后, 厚度的变化不再对计算得出的尺度因子产生明显的影响; 另一方面, 说明我们的假设并未完全考虑所有可能对有效基质扩散系数产生影响的因素, 可能存在其他物理机制造成的有效基质扩散系数的尺度效应, 如岩石基质的非均质性、裂缝表面的粗糙程度、统计结果的精确性以及达西定律的适用性等因素。
图7 尺度因子与观测尺度的关系
Fig. 7 Relationship between the scale factor and observation scale
本文通过建立考虑裂缝周围损伤区的三维数值模型, 模拟裂缝及基质中的溶质运移过程, 分析裂缝周围损伤区对基质扩散的影响, 提出一种解释有效基质扩散系数尺度效应的机理。主要结论如下。
1)裂缝周围损伤区的厚度、孔隙度及渗透率均会对有效基质扩散系数产生影响。其中, 随着损伤区厚度及孔隙度的增加, 有更多的溶质进入岩石基质并在其中停留, 增强了岩石中的分子扩散过程, 使得测量得到的有效基质扩散系数增大。损伤区厚度对有效基质扩散系数的影响在 0.01～0.1m 之间, 存在一个上限。
2)大尺度裂缝或断层周围的损伤区存在更密集的小尺度裂缝, 该区域岩石比远离裂缝面的未破坏岩石具有更大的孔隙度和渗透率。随着裂缝尺度增加, 损伤区厚度、渗透率及孔隙度均会增加, 从而增大有效基质扩散系数。
3)有效基质扩散系数的尺度因子与观测尺度之间存在正向关系, 其斜率为 0.29, 与现场试验数据之间存在一定的差距。因此, 裂缝周围存在损伤区的假设只能部分地解释有效基质扩散系数的尺度效应。
由于地下裂隙介质的复杂性, 可能有多种因素造成有效基质扩散系数的增大。如岩石基质可能具有不同程度的非均质性, 地下裂缝表面通常具有一定的粗糙度, 这些因素都可能造成溶质运移时间的增加, 使得计算得出的有效基质扩散系数增加。未来的研究中需要考虑这些因素对有效基质扩散系数的影响。
参考文献
[1] Bodin J, Delay F, de Marsily G. Solute transport in a single fracture with negligible matrix permeability: 1. Fundamental mechanisms. Hydrogeology Journal, 2003, 11(4): 418–433
[2] Maloszewski P, Zuber A. Mathematical modeling of tracer behavior in short-term experiments in fissured rocks. Water Resources Research, 1990, 26(7): 1517–1528
[3] Callahan T J, Reimus P W, Bowman S, et al. Using multiple experimental methods to determine fracture/ matrix interactions and dispersion of nonreactive solu-tes in saturated volcanic tuff. Water Resources Resear-ch, 2000, 36(12): 3547–3558
[4] Reimus P W, Callahan T J. Matrix diffusion rates in fractured volcanic rocks at the Nevada Test Site: evidence for a dominant influence of effective frac- ture apertures. Water Resources Research, 2007, 43: W07421
[5] Zhou Q, Liu H, Molz F J, et al. Field-scale effective matrix diffusion coefficient for fractured rock: results from literature survey. Journal of Contaminant Hydro-logy, 2007, 93: 161–187
[6] Shapiro A M. Effective matrix diffusion in kilometer-scale transport in fractured crystalline rock. Water Re-sources Research, 2001, 37(3): 507–522
[7] Neretnieks I. A stochastic multi-channel model for so-lute transport-analysis of tracer tests in fractured rock. Journal of Contaminant Hydrology, 2002, 55: 175–211
[8] Dai Z, Wolfsberg A, Lu Z, et al. Upscaling matrix di-ffusion coefficients for heterogeneous fractured rocks. Geophysical Research Letters, 2007, 34(7): L07408
[9] Lei D, Sun H, Zhang Y, et al. Upscaling solute transport in rough single-fractured media with matrix diffusion using a time fractional advection-dispersion equation. Journal of Hydrology, 2023, 627: 130280
[10] Liu H, Zhang Y Q, Zhou Q, et al. An interpretation of potential scale dependence of the effective matrix diffusion coefficient. Journal of Contaminant Hydrolo-gy, 2007, 90(1/2): 41–57
[11] Evans J P, Forster C B, Goddard J V. Permeability of fault-related rocks, and implications for hydraulic structure of fault zones. Journal of Structural Geolo-gy, 1997, 19(11): 1393–1404
[12] Walker R J, Holdsworth R E, Armitage P J, et al. Fault zone permeability structure evolution in basalts. Geo-logy, 2013, 41(1): 59–62
[13] Choi J H, Edwards P, Ko K, et al. Definition and classification of fault damage zones: a review and a new methodological approach. Earth-Science Reviews, 2016, 152: 70–87
[14] Mitchell T M, Faulkner D R. Towards quantifying the matrix permeability of fault damage zones in low poro-sity rocks. Earth and Planetary Science Letters, 2012, 339: 24–31
[15] Webster D S, Proctor J F, Marine I W. Two-well tracer test in fractured crystalline rock (No. 1544-I). Wa-shington: USGPO, 1970
[16] Tang D H, Frind E O, Sudicky E A. Contaminant tran-sport in fractured porous media: analytical solution for a single fracture. Water Resources Research, 1981, 17 (3): 555–564
[17] Liou T. Interpretation of the enhancement of field-scale effective matrix diffusion coefficient in a single frac-ture using a semi-analytical. Hydrogeological Proces-ses, 2009, 23: 816–829
[18] Lamur A, Kendrick J E, Eggertsson G H, et al. The permeability of fractured rocks in pressurised volcanic and geothermal systems. Scientific Reports, 2017, 7 (1): 6173
[19] Torabi A, Berg S S. Scaling of fault attributes: a review. Marine and Petroleum Geology, 2011, 28(8): 1444–1460
[20] Bense V F, Gleeson T, Loveless S E, et al. Fault zone hydrogeology. Earth-Science Reviews, 2013, 127: 171–192
Effect of Fracture Damage Zones on the Effective Matrix Diffusion Coefficient of Fractured media
Abstract Field tracer tests indicate that the effective matrix diffusion coefficient of fractured media increases with the increase of observation scale. To understand the mechanisms of this phenomenon, a 3D fracture-matrix model considering a horizontal fracture as well as the damage zone around the fracture is developed. The effect of the damage zone properties, including width, porosity and permeability, on the effective matrix diffusion coefficient is analyzed through a series of solute transport simulations. The results show that when the damage zone has a larger thickness and porosity, more solute enters the rock matrix, which leads to an increase in the measured effective matrix diffusion coefficient. For natural fractures, the permeability and porosity of damage zone increase with the decrease of the distance from the fracture surface. With the increase of fracture scale, the thickness, permeability and porosity of the damage zone increase. The simulation results of different observation scales show that a positive relationship exists between the scale factor of the effective matrix diffusion coefficient and observation scale. This study interprets part of the mechanism of the scale dependence of the coefficient and contributes to the understanding of the solute transport process in complex fracture-matrix systems.
Key words matrix diffusion; scale dependence; damage zone; solute transport; fractured media