Reactive Transport Numerical Modeling of Typical Heavy Metal Pollutants in Three-Dimensional Fracture Networks
-
摘要: 自然界中裂隙介质的具体结构往往控制着其中的水流流动和溶质运移过程. 由于真实裂隙网络的复杂性，常据不同裂隙参数（形状、位置、大小、密度、走向等）概化并使用特定数学分布来刻画和产生离散裂隙网络，但这种离散裂隙网络介质因其高度非均质性和各向异性导致传统基于多孔介质理论的数值方法难于模拟裂隙水流和溶质运移过程. 利用三维离散裂隙网络模型生成场地裂隙介质不连续体，并基于裂隙空间参数将离散裂隙网络模型映射为等效多孔介质模型，由此进一步耦合地下水流与多组分反应性溶质运移开源代码PFLOTRAN，以实现三维裂隙网络中重金属污染物反应运移过程的数值模拟与污染物通量的定量评估. 算例研究表明，随着裂隙密度的增大，裂隙网络表征逐渐趋近于多孔介质；裂隙发育的密度和方向均影响着裂隙网络的连通性，进而影响裂隙水流中不同组分的迁移规律和通量评估，其中保守性组分与不同反应性组分的迁移行为有明显差异. 同时，研究表明，裂隙粗糙度对刻画裂隙网络中的污染组分迁移过程和通量评估具有重要影响. 本文方法可为场地裂隙介质中地下水流和污染物的反应运移提供有效模拟手段和定量评估工具.Abstract: The geologic discontinuities present in natural fractured rock governs groundwater flow and solute transport process. Due to the complexity of real-world fracture network, discrete fracture network (DFN) is often be described and generated based on specific probability distributions determined by different parameters (shape, location, size, density, orientation, etc.). Due to DFN's highly heterogeneous and an isotropic features, it's a challenge to simulate groundwater flow and solute transport processes using traditional numerical solution based on the experience achieved in porous media. In this paper, a 3-D DFN model is built for field-scale study, and mapped to an equivalent porous medium (EPM) model, which further couples groundwater flow model with the multi-component reactive transport code PFLOTRAN to achieve numerical simulation of the reactive transport process of heavy metal contaminants in the fractured network and quantitative assessment of mass fluxes. The study shows that the behavior of the fracture network gradually transforms into the transport in porous media as the density of the network increases; Both the density and orientation of fracture development affect the connectivity of the fracture network, which in turn affects the migration pattern and flux assessment of different components in fractured rock groundwater flow, where the migration behavior of conservative components and various reactive components are significantly different.Meanwhile, it is shown that roughness has an important influence on the migration process and flux assessment of the transport components in the fracture network.This paper provides an effective simulation and quantitative assessment tool for the reactive transport of heavy metals and groundwater flow in the fractured field.
-
0. 引言
水文地质学中将节理、局部断裂等面状不连续结构统一称之为裂隙（陈宏峰等，2016）. 自然界中裂隙尺度可从厘米级到十千米级，为研究裂隙中的溶质运移，首先需准确刻画裂隙网络，并掌握其中流体的运动规律. 这是其区别于多孔介质中溶质运移的研究难点.
从单一裂隙到裂隙网络，由理想平行板模型到真实曲折裂隙空间，是裂隙研究的理论渐进过程. 裂隙中地下水流动速度一般较慢，Reynolds数较小，因此可建立基于达西定律的水流运动偏微分方程并由此数值求解. 这种简化涉及到实际岩体的典型单元体表征，得到的是平均化的参数；广泛使用的立方定律是基于光滑平行板假设并通过简化Navier-Stokes方程推导得到的（Zimmerman and Yeo, 2000），而自然界中的裂隙并不一定符合该假设（李一鸣和文章，2020）. 通过在小尺度上直接求解Navier-Stokes方程来修正立方定律（Trinchero et al.，2021）或采用实验方法拟合非达西流动的具体参数（甘磊等，2021）是当今多数研究的主要手段.
另一方面，现实中如坝基渗流、水力压裂、核废料储存等实际应用均涉及到裂隙与应力场之间的关系以及渗流过程（王伟等，2021），这又促使了理论在场地尺度的应用推广，包括场地试验验证与模型参数率定（Wilson et al.，1983）. 例如，Long et al.（1990）在Stripa场地使用场地实测数据建立了包含主要裂隙分布特征的水文地质概念模型；Peters and Klavetter，（1988）在Yucca使用连续介质模型评估了裂隙岩体中的非饱和渗流运动. 使用离散裂隙网络（Discrete Fracture Network，DFN）以进行地质力学表征及水力学行为考察的具体数值研究已有多种不同的实现（有限元、离散元等），基于多物理过程的地质力学模拟耦合裂隙岩体中水力学过程以生成更具有天然裂隙网络的建模流程也在不断发展中（Lei et al.，2017）. 在获取了生成的裂隙网络后，由于裂隙与岩体骨架之间的巨大差异，学者们（Li et al.，2020；Hu et al.，2022）通过引入裂隙岩体和骨架中的交换项以考察三维裂隙网络及其骨架中溶质运移问题，但一般空间范围（长宽高）不超过十几米.
在较小尺度上，可对岩芯样品直接采用断层扫描或其他技术获取真实的裂隙分布（Yang et al.，2021），而在较大尺度上，可利用地震波、野外踏勘，结合卫星图片解译等方法概略厘定大地构造和宏观裂隙的基本发育情况，此时裂隙分布是确定性的. 但在场地尺度上，利用测井、钻孔、测窗等测量手段获取完整精细的裂隙网络是非常困难的，为此采用DFN模型来刻画场地尺度的裂隙分布是一种可行的途径.
由于DFN只考虑裂隙内的地下水流动，在模拟溶质运移时不能反应岩石骨架中的水动力弥散机制. 因此，本文采用将DFN映射至等效多孔介质（rquivalent porous media，EPM）的方法来进行计算. 在映射过程中，实际上是对裂隙网络进行了一定程度的升尺度处理，并使用传统的对流-弥散方程描述溶质运移过程，有效降低了计算量. 在保证裂隙网络中地下水流场模拟精度的前提下进行地下水溶质的反应性迁移模拟，为类似裂隙发育的场地地下水污染过程预测评价提供了一种高效的数值模拟方法.
1. 研究方法
1.1 离散裂隙网络模型(DFN)
描述裂隙的具体参数有走向、倾角、大小、位置和形状等（Makedonska et al.，2016）. 本文采用提供椭圆盘模型（Dershowitz and Einstein, 1988）的DFNWORKS来生成DFN，通过指定一系列二维平行板构成三维裂隙网络（Hyman et al.，2014）. 其使用一种特征抑制网格划分算法（feature rejection algorithm for meshing，FRAM），以生成逼近用户设定的目标裂隙数量或P32值（裂隙面积与生成空间体积的比值，单位为m-1）的裂隙网络，并保证网络中的所有裂隙特征都在某一截断长度之上，并可以将与生成裂隙网络空间（六面体）中指定面不连通的裂隙面移除. 本文中使用DFNWORKS生成的DFN包括以下假设：
（1）生成的单个裂隙形状为椭圆形（当椭圆的长短轴之比为1时则为圆形）.
（2）生成裂隙的形心在三维空间中随机均匀分布，并可移除不连通的裂隙.
（3）生成裂隙的朝向服从von Mises-Fisher分布，其分布函数满足：
f(α)=κ⋅sinα⋅eκ⋅cosαeκ−e−κ, (1) 式中：α为裂隙极点与平均方位的偏差；参数κ为表征聚集程度. 自然界中不同场地κ值有所不同，在某些场地根据测线上观测裂隙的赤平投影图统计，通常κ值在10~20（Romano et al.，2020）. 在生成裂隙时，本文取κ=10来反映裂隙的集中程度.
（4）生成裂隙的大小（半径）服从截断幂律分布，其采样半径R可表示为：
R=R0⋅[1−u+u⋅(R0Ru)γ]−1γ， (2) 式中：R0和Ru分别为截断值的下界和上界（m）；u为在0到1间的随机数；γ为衰减因子，可取γ=2（Makedonska et al.，2016）.
（5）生成裂隙的隙宽与裂隙大小相关，可表示为：
b=F⋅Rβ， (3) 式中：F和β为相关关系参数，本文分别取值5×10-5和0.5（Hyman et al.，2016）.
DFNWORKS只对裂隙网络本身使用三角网进行空间离散（Hyman et al.，2015），并可采用粒子追踪的方法进行溶质运移计算（Makedonska et al.，2015）. DFNWORKS还可利用八叉树加密网格剖分算法（Sweeney et al.，2020），由加密的四面体完全非结构化网格来考虑岩石骨架的影响. 由于相邻网格中心和公用面之间的非正交性，使用两点流格式（two-point flux-approximation）的数值计算代码将引入计算误差，加密等级较高时，网格数量导致计算负荷也较大.
1.2 映射至等效多孔介质模型(EPM)
给定DFN中每个裂隙都为椭圆形的二维平行板. 若裂隙的宽度为$ {b}_{f} $，其映射过程就是将渗透系数张量进行旋转到[0，0，1]方向. 对于速度矢量$ \vec{v} $以及水力梯度矢量$ \vec{J} $，在旋转后的新坐标系记为$ {\vec{v}}^{\prime} $、$ {\vec{J}}^{\prime} $：
{→v′=R·→v→J′=R·→J. (4) 在每个裂隙中，旋转前后的介质中有达西定律：
{→v=[K]→J→v′=[K]→J′. (5) 若裂隙平面法向向量（单位向量）为$ \vec{n} $，记[0，0，1]为$ \vec{a} $，则可求得$ \vec{n} $与$ \vec{a} $之间的旋转夹角θ：
θ=arccos(→a·→n|a|⋅|n|). (6) 令$ \vec{k}=\vec{a}\times \vec{n} $，该叉乘结果$ \vec{k} $即为旋转轴. 由罗德里格斯旋转公式，并令$ {w}_{1}=\mathrm{c}\mathrm{o}\mathrm{s}\frac{\theta }{2} $，$ {w}_{2}=\mathrm{s}\mathrm{i}\mathrm{n}\frac{\theta }{2}\cdot {k}_{x} $，$ {w}_{3}=\mathrm{s}\mathrm{i}\mathrm{n}\frac{\theta }{2} $，$ {w}_{4}=\mathrm{s}\mathrm{i}\mathrm{n}\frac{\theta }{2}\cdot {k}_{z} $，则旋转矩阵R为：
[1−2w23−2w242w2w3−2w1w42w2w4+2w1w32w2w3+2w1w41−2w22−2w242w3w4−2w1w22w2w4−2w1w32w3w4+2w1w21−2w22−2w23]. (7) 计算R·RT，根据t2+u2+v2+w2=1化简可算出结果为三阶单位矩阵，可知旋转矩阵R为正交矩阵，有R-1=RT. 将式（4）~（7）及正交矩阵特性代入整理得到映射前后渗透率张量的关系：
[K]′=R⋅[K]⋅RT. (8) 旋转前单个裂隙渗透张量根据立方定律计算，即$ {k}_{f}=\frac{{b}_{f}^{2}}{12} $，旋转后计算其固有导水系数Tf （m3），kf'为根据式（10）旋转后的裂隙渗透张量， Tf=bfkf'. 将Tf除以结构化网格边长d，即完成升尺度映射，只取主对角线的元素：
[kxx000kyy000kzz]=1d[Txx000Tyy000Tzz]. (9) 最后将空间中各个裂隙的贡献累加，岩石骨架渗透率统一取为10-20 m2：
k=∑Mm=1kmf, (10) 式中：$ {k}_{f}^{m} $表示网格中第m条裂隙的渗透率；M是网格中裂隙的总数.
这种映射方法将原本作为“零厚度”的裂隙网络提升到了网格尺度. 图 1展示了EPM中的映射裂隙部分，在未映射前得到的裂隙网络为使用三角单元组成的非结构化网格，映射后变为每个单元大小均等的结构化网格. 利用开源脚本Mapdfn可实现以上DFN映射到EPM的计算过程，进而获得对应裂隙岩体的渗透系数张量（钟启明等，2006）.
1.3 三维裂隙水流与溶质反应运移模型（PFLOTRAN）
在PFLOTRAN中单相变饱和等温系统的控制方程如下，模拟中水相的饱和度为1，温度为25℃.
∂(φρ)∂t+∇⋅(ρq)=0, (11) q=−kμ∇(P−ρgz), (12) $ \mathrm{式}\mathrm{中}:\varphi $为孔隙度（无量纲）；$ \rho $为流体密度（kg$ \cdot $m-3）；q为达西流速（m$ \cdot $s-1）；k为渗透率（m2）；$ \mu $为流体动力粘滞系数（Pa$ \cdot $s）；P为压力（Pa）；g为重力加速度（m$ \cdot $s-2）；z为位置高程（m）；裂隙介质的孔隙度可通过累加穿过每个网格中裂隙的隙宽贡献来计算（岩石骨架孔隙度可统一取值0.005）：
φ=1d∑Mm=1bmf, (13) 式中：$ {b}_{f}^{m} $表示网格中第m条裂隙的隙宽，其他符号同上.
对于多组分溶质运移模型，应用质量守恒定律可得到组分j的对流-弥散方程：
∂(φΨj)∂t+∇⋅(q−φD∇)Ψj=Qj−∑mυjmIm, (14) 式中：$ {\varPsi }_{j} $表示组分j以自由离子和络合离子存在的总浓度（mol$ \cdot $m-3）；$ {Q}_{j} $为源汇项（mol$ \cdot $m-3$ \cdot $s-1）；$ {I}_{\mathrm{m}} $为矿物反应速率（mol$ \cdot $m-3$ \cdot $s-1）；$ {\upsilon }_{j\mathrm{m}} $表示对组分j有贡献的化学计量数（无量纲）；D为水动力弥散张量（m2$ \cdot $s-1）. 在PFLOTRAN中，计算D时只取其主对角线（Pandey and Rajaram, 2016），即：
[τDm+αLv2x|v|+αTHv2y|v|+αTVv2z|v|000τDm+αLv2y|v|+αTHv2x|v|+αTVv2z|v|000τDm+αLv2z|v|+αTVv2x|v|+αTVv2y|v|], (15) 式中：$ {D}_{m} $为分子扩散系数（m2$ \cdot $s-1）；τ为曲迂度，与孔隙度相关，本文中表示为τ=λ/φ，λ为人为指定参数，最终τ的范围在0.2~3；$ {\alpha }_{\mathrm{L}} $为纵向弥散度（m）；$ {\alpha }_{\mathrm{T}\mathrm{H}} $为水平横向弥散度（m）；$ {\alpha }_{\mathrm{T}\mathrm{V}} $为垂直横向弥散度（m）. 本文中取$ {\alpha }_{\mathrm{L}} $=1，$ {\alpha }_{\mathrm{T}\mathrm{H}}=0.005 $、$ {\alpha }_{\mathrm{T}\mathrm{V}}=0.005 $；$ {v}_{x} $、$ {v}_{y} $、$ {v}_{z} $分别为速度矢量在x、y、z方向的分量，$ \left|v\right| $为模长.
对于反应性矿物组分来说，其反应性动力学速率为：
∂φm∂t=¯VmIm, (16) 式中：$ {\varphi }_{\mathrm{m}} $为矿物体积分数（无量纲）；$ \overline{{V}_{\mathrm{m}}} $为矿物的摩尔体积（m3$ \cdot $mol-1）；其中总浓度$ {\varPsi }_{j} $可进一步表示为：
Ψj=Cj+∑Nxjk=1υjiCi, (17) 式中：$ {C}_{j} $组分j的自由离子浓度（mol$ \cdot $m-3）；$ {C}_{i} $与组分j有关的次级离子i的浓度（mol$ \cdot $m-3）；$ {\upsilon }_{ji} $为其次级反应的化学计量数（无量纲）；$ {N}_{xj} $络合离子的种类个数（无量纲）；
次级组分i的浓度$ {C}_{i} $可根据平衡反应的质量作用定律计算得到：
Ci=Kiγi∏j(γjCj)υji, (18) 式中：$ {K}_{i} $为反应平衡常数；$ {\gamma }_{i}\mathrm{为} $活度系数（无量纲）.
矿物反应动力学速率可表示为（Lichtner，2016）：
Im=−kmAm(1−QmKm), (19) $ \mathrm{式}\mathrm{中}:{k}_{\mathrm{m}} $为反应速率常数（mol$ \cdot $m-2$ \cdot $s-1）；$ {A}_{\mathrm{m}} $为矿物比表面积（m-1）；$ {K}_{\mathrm{m}} $为溶度积；$ {Q}_{\mathrm{m}} $为活度积.
2. 研究算例分析
2.1 场地背景与概念模型
本文方法拟应用于一个酸性尾矿向三维裂隙介质中释放污染物的运移过程模拟（Molson et al.，2012）. 如图 2所示，模型空间概化为200 m×100 m×50 m的长方体区域，地表有1.25 m的均匀盖层（图示黄色区域），其孔隙度为0.1，曲迂度为0.2，水平渗透率为10-11 m2，垂直渗透率为10-12 m2. 模拟区域为完全饱和流动，地下水初始平均水力梯度为0.002，西高东低. 其中西部边界为水头50 m的定水头边界，东部为通用水头边界，南边边界为零通量边界，底部为隔水边界；区域顶部接受恒定的面状补给，其补给强度为200 mm/a. 图 2中青色区域为下文中穿透曲线及通量计算的考察断面（不包括盖层），计算通量时取断面上的浓度和速度的平均值.
假设组成裂隙岩石的矿物有7种，包括方解石、菱铁矿、三水铝石、水铁矿、软石膏、铅矾和白铅矿，溶质运移计算中涉及13种主要离子（H+、Al3+、Fe3+、Fe2+、Mn2+、SO42-、CO32-、Cl-、K+、Na+、Mg2+、Ca2+、Pb2+参与对流弥散方程运算）和3种次级离子（OH-、HCO3-、CO2（aq））. 污染源（图中表层红色区域）空间所在位置范围：x=90~110 m，y=48~52 m，z=48.75~49.00 m. 初始浓度、面状补给中各组分浓度及尾矿区域源强均为给定条件，具体取值据文献，参见表 1（Walter et al.，1994）. 为削弱数值弥散带来的影响，减小Peclet数，模型空间离散剖分的网格大小为0.625 m×0.625 m×0.625 m.
表 1 溶质运移模型初始及边界各组分浓度（单位：mol$ \cdot $L-1, pH除外）Table Supplementary Table Initial and boundary condition of the solute transport model (unit: mol$ \cdot $L-1, except pH)主要离子组分 背景值 源强组分 降水组分 H+ (pH) 9.2 4.0 7.0 Al3+ 1.00×10-20 4.30×10-3 1.30×10-7 Fe3+ 1.00×10-22 2.00×10-7 2.30×10-8 Fe2+ 4.20×10-5 3.06×10-2 5.40×10-5 Mn2+ 2.90×10-5 7.84×10-3 4.70×10-5 SO42- 5.30×10-3 5.00×10-2 7.50×10-3 CO32- 2.10×10-5 4.92×10-4 3.90×10-3 Cl- 1.00×10-3 1.58×10-2 1.00×10-3 K+ 6.60×10-5 8.14×10-4 6.70×10-5 Na+ 1.30×10-3 1.39×10-3 1.30×10-3 Mg2+ 1.50×10-3 9.69×10-4 1.90×10-3 Ca2+ 5.20×10-3 1.08×10-2 6.09×10-3 Pb2+ 1.00×10-8 1.52×10-5 1.00×10-20 2.2 裂隙网络参数与情景设置
对于给定分布的DFN而言，只有生成裂隙达到一定的数量（密度），才能联通区域中不同的裂隙面从而形成裂隙网络. 一般可使用逾渗阈值（percolation parameter）（de Dreuzy et al.，2012）来衡量这种全局连通性. 在裂隙服从分布及生成空间范围确定的情况下，可直接给出逾渗阈值的计算公式，也可通过简单的试错法来估算产生连通裂隙网络的裂隙数量大致范围. 基于本文场地范围，当目标裂隙数量在250左右时，可生成连通的裂隙网络，而当分别取2倍和4倍于该数值的裂隙数量时则可生成更为密集的裂隙网络，对应这3种不同情境的裂隙网络分别标记为低密度、中密度、高密度裂隙网络. 不同裂隙网络的相关参数和移除相应不连通裂隙的网络可分别参见表 2和图 3所示. 所有裂隙网络的倾角均设定为90°.
表 2 不同密度连通裂隙网络的参数Table Supplementary Table Parameters of connected fracture networks with different density裂隙密度 目标裂隙数N/(条) 连通裂隙数n/(条) P32 低密度 250 54 0.16 中密度 500 277 0.70 高密度 1 000 511 1.12 对于不同的裂隙密度，随着裂隙方向的不同，生成的裂隙网络亦随之不同. 当中密度（N=500）裂隙网络分别将其走向90°改变为60°、75°、105°和120°时，生成的连通裂隙网络相关参数如表 3所示.
表 3 中密度（N=500）条件下生成不同走向的裂隙网络特征Table Supplementary Table Fracture network features with different trends generated under the condition of medium density (N=500)不同情景 走向φ(°) 连通裂隙数n(条) P32 Case 0（基准方向） 90 277 0.70 Case 1 60 192 0.51 Case 2 75 92 0.24 Case 3 105 114 0.32 Case 4 120 256 0.66 另一方面，裂隙的粗糙度影响其渗透性，因此许多学者通过各种单裂隙试验来修正立方定律（Barton et al.，1985），本文采用Louis对立方定律的修正公式（梁敏，2010）：
T=b312⋅11+8.8(Δ2b)1.5, (20) 式中：∆为单条裂隙中凸起的绝对高度；b为最大裂隙宽度.
对于高密度（N=1 000）的裂隙网络，公式（20）中$ \frac{\Delta }{b} $分别取值为1、0.5、0（未修正的立方定律）时，则可得到渗透率不同但几何拓扑关系相同的3个裂隙网络，分别记为低渗透率、中渗透率、高渗透率裂隙网络. 它们之间的渗透率是一个简单的倍数关系，映射后对应裂隙网络的kxx、kyy、kzz的范围参见表 4.
表 4 不同渗透率裂隙网络设定Table Supplementary Table Fracture networks settings with different permeability渗透特性 修正系数 渗透率取值范围(m2) $ \overline{{\boldsymbol{k}}_{\boldsymbol{x}\boldsymbol{x}}} $ $ \overline{{\boldsymbol{k}}_{\boldsymbol{y}\boldsymbol{y}}} $ $ \overline{{\boldsymbol{k}}_{\boldsymbol{z}\boldsymbol{z}}} $ 低渗透率 0.24 [2.95×10-13, 3.58×10-12] [5.32×10-17, 1.54×10-12] [3.18×10-13, 3.47×10-12] 中渗透率 0.48 [5.78×10-13, 7.02×10-12] [1.04×10-16, 3.02×10-12] [6.22×10-13, 6.79×10-12] 高渗透率 1.00 [1.21×10-12, 1.47×10-11] [2.19×10-16, 6.35×10-12] [1.31×10-12, 1.43×10-11] 2.3 重金属组分与反应运移参数
前已述及，组成裂隙岩石的矿物有方解石、菱铁矿、三水铝石、水铁矿、软石膏、铅矾和白铅矿等7种，主要的动力学溶解沉淀反应参见表 5.
表 5 化学反应及其平衡常数Table Supplementary Table Chemical reactions and their equilibrium constants反应表达式 平衡常数lgK $ \mathrm{O}{\mathrm{H}}^{-}+{\mathrm{H}}^{+}\leftrightarrow {\mathrm{H}}_{2}\mathrm{O} $ 13.995 1 $ {\mathrm{H}}^{+}+\mathrm{C}{\mathrm{O}}_{3}^{2-}\leftrightarrow \mathrm{H}\mathrm{C}{\mathrm{O}}_{3}^{-} $ 10.328 8 $ \mathrm{C}{\mathrm{O}}_{2}\left(\mathrm{a}\mathrm{q}\right)+{\mathrm{H}}_{2}\mathrm{O}\leftrightarrow {\mathrm{H}}^{+}+\mathrm{H}\mathrm{C}{\mathrm{O}}_{3}^{-} $ -6.344 7 $ \mathrm{C}\mathrm{a}\mathrm{C}{\mathrm{O}}_{3}+{\mathrm{H}}^{+}\leftrightarrow \mathrm{C}{\mathrm{a}}^{2+}+\mathrm{H}\mathrm{C}{\mathrm{O}}_{3}^{-} $ 1.848 7 $ \mathrm{F}\mathrm{e}\mathrm{C}{\mathrm{O}}_{3}+{\mathrm{H}}^{+}\leftrightarrow \mathrm{F}{\mathrm{e}}^{2+}+\mathrm{H}\mathrm{C}{\mathrm{O}}_{3}^{-} $ -0.192 0 $ \mathrm{A}\mathrm{l}{\left(\mathrm{O}\mathrm{H}\right)}_{3}+3{\mathrm{H}}^{+}\leftrightarrow \mathrm{A}{\mathrm{l}}^{3+}+3{\mathrm{H}}_{2}\mathrm{O} $ 7.756 0 $ \mathrm{F}\mathrm{e}{\left(\mathrm{O}\mathrm{H}\right)}_{3}+3{\mathrm{H}}^{+}\leftrightarrow \mathrm{F}{\mathrm{e}}^{3+}+3{\mathrm{H}}_{2}\mathrm{O} $ 4.896 0 $ \mathrm{C}\mathrm{a}\mathrm{S}{\mathrm{O}}_{4}\cdot 2{\mathrm{H}}_{2}\mathrm{O}\leftrightarrow \mathrm{C}{\mathrm{a}}^{2+}+\mathrm{S}{\mathrm{O}}_{4}^{2-}+2{\mathrm{H}}_{2}\mathrm{O} $ -4.482 3 $ \mathrm{P}\mathrm{b}\mathrm{S}{\mathrm{O}}_{4}\leftrightarrow \mathrm{P}{\mathrm{b}}^{2+}+\mathrm{S}{\mathrm{O}}_{4}^{2-} $ -7.852 7 $ \mathrm{P}\mathrm{b}\mathrm{C}{\mathrm{O}}_{3}+{\mathrm{H}}^{+}\leftrightarrow \mathrm{P}{\mathrm{b}}^{2+}+\mathrm{H}\mathrm{C}{\mathrm{O}}_{3}^{-} $ -3.209 1 表 5中化学反应的平衡常数均来源于PFLOTRAN数据库（Hammond and Lichtner, 2010）. 岩石骨架中方解石、菱铁矿、三水铝石、水铁矿的初始体积分数为0.001，其他矿物的初始体积分数均为0. 方程式（19）中，方解石的反应速率常数（km）为10-8 mol$ \cdot $m-2$ \cdot $s-1，其他各矿物的反应速率常数（km）均为10-10 mol$ \cdot $m-2$ \cdot $s-1（Molson et al.，2008），据实验室测定（Jordan and Rammensee, 1998）并减小以适应于场地尺度（White and Brantley, 2003）矿物的比表面积均给定为1 m-1. 模型运行在25 ℃等温条件下，模拟时长为1 000 d.
3. 结果与讨论
模拟结果显示，相同的组分在不同裂隙网络中具有类似的迁移和分布规律. 以高密度（高渗透率）裂隙网络1 000 d时的模拟结果为例，对于作为保守性溶质（如Cl-）来说，主要在裂隙网络中的优势通道中运移（图 4a中对应高流速区），并沿岩体骨架有一定范围的弥散（如图 4a中截面中的浓度分布）；但反应性组分（如Fe2+）的运移过程与保守性组分相比有着显著的不同（如图 4b）. 为便于比较，图 4b中所示的Cl-和菱铁矿（FeCO3）饱和指数分布分别向两侧平移35 m（实际上Cl-、Fe2+以及FeCO3的饱和指数大致集中在同一个空间范围，为方便观察平移），并在菱铁矿饱和指数分布中用黄色绘制pH=6.5的锋面，由于菱铁矿的沉淀溶解反应主要受pH控制，导致反应性Fe2+较保守性Cl-的迁移速度更为缓慢. 从菱铁矿饱和指数的空间分布可看出，在酸性峰面以外是碳酸亚铁主要沉淀（生成菱铁矿）的区域，而在pH较低区域内（污染源区），FeCO3饱和指数小，菱铁矿（包括方解石）都溶解，三水铝石和水铁矿会形成微量的沉淀.
图 4 高密度裂隙网络中不同组分运移结果及其空间分布a.对流速度场分布及保守性Cl-浓度空间分布；b. Cl-、Fe2+浓度以及菱铁矿饱和指数空间分布；c.Pb2+浓度及碳酸铅、硫酸铅、软石膏饱和指数的空间分布Fig. 4. Solute transport results and their spatial distributions of different components in the high-density fracture network对于反应性组分Pb2+来说，它在运移过程中可同时生成PbSO4与PbCO3两种沉淀（如图 4c所示，类似图 4b中将结果向左右平移以方便展示）. 在该体系中，Ca2+和SO42-在远离源区一定距离后，两者在合适条件下会部分沉淀成为软石膏（CaSO4·2H2O），同时PbCO3的沉淀也降低了远离源区的Pb2+浓度，这些直接反应的竞争间接地导致PbSO4只能在接近源区产生沉淀. 而PbCO3发生沉淀的空间位置则与FeCO3类似. 由此可见，反应性重金属组分Pb2+的迁移过程是多种水文地球化学反应共同作用的结果.
3.1 裂隙密度对运移结果的影响
如表 2所示，当保持裂隙方向不变而裂隙密度增大时，连通的裂隙数量增多（裂隙网络的连通性增强），介质的非均质性相对减小，且当裂隙密度增大到一定程度时，其总体裂隙网络趋近于多孔介质. 当裂隙网络密度过高时，由于离散裂隙网络映射到等效多孔介质模型算法的特点，决定了低渗透性的岩石骨架将基本缺失，除非将结构化网格的剖分网格足够小，但由此会带来巨量网格而失去这种映射方法的优势，此时本文映射方法是不合适的.
随着裂隙密度的增高，由于裂隙网络的连通性变好，对于保守性Cl-来说，其在下游断面上的平均浓度（图 5a）和最大浓度均随之增大，同时其通量亦变大；但反应性组分来说，随着裂隙网络密度的增加，最终下游断面的平均浓度虽然总体呈现为增加的态势，但其最大浓度则表现不同，即由于不同裂隙网络对化学反应的影响，使得通过下游断面的反应性组分的最大浓度并不总是随着裂隙数量的增加而单调减小，而是表现为上下浮动. 如图 5b所示，反应性Fe2+在高密度情景下通过下游断面的最大浓度小于低密度情景下的最大浓度，说明在低密度的裂隙网络中，优势流动的现象更为明显（Tsang and Tsang, 1987）. 当裂隙网络密度低时，溶质运移集中于几个大的主干裂隙中，从而使得下游断面的最大浓度更高，与此同时，主干中的高浓度溶质又会对沉淀溶解造成影响，从而又降低了反应性溶质的浓度，这也是亚铁离子最终最大浓度在中密度和高密度网络中相差不大的原因.
图 5 不同密度裂隙网络情景下不同组分在下游断面的浓度穿透曲线a. Cl-平均浓度；b. Fe2+最大浓度Fig. 5. Breakthrough curves of different components in the downstream section under different density fracture network scenarios3.2 裂隙方向对运移结果的影响
裂隙网络中溶质的迁移结果不仅受到裂隙密度的影响，还受控于裂隙发育的方向. 在同一裂隙密度条件下，不同方向的裂隙介质形成的裂隙网络具有不同的连通性. 如表 3所示，当改变中密度裂隙网络的裂隙走向时，所生成的裂隙网络中对应的连通裂隙数量不同，φ=90°方向对应的连通裂隙最多，P32值最大；60°和120°方向对应的连通裂隙数与其P32值次之；而75°和105°方向对应的连通裂隙数与其P32值最小. 也就是说，随着裂隙方向的改变，生成裂隙网络的连通性跟着变化，进而导致下游断面溶质的对应穿透曲线亦随之变化. 如图 6a，6b所示，无论是Cl-还是Fe2+，都是90°方向对应裂隙网络出流断面的溶质穿透曲线在1 000 d时的浓度最高，75°方向情景的浓度最低，其他方向同样随着连通性的增强其浓度随之升高. 尤其对于90°的裂隙网络来说，溶质组分由源到汇（出流断面）的流程最短，故而其浓度最高.
图 6 不同走向裂隙网络情景下不同组分在下游断面的穿透曲线a. Cl-平均浓度；b. Fe2+平均浓度；c. Pb2+平均浓度；d. Pb2+通量Fig. 6. Breakthrough curves of different components in the downstream section under different trends fracture network scenarios然而，对于Pb2+来说，如图 6c所示，φ=120°裂隙网络下游断面的浓度反而高于φ=90°的裂隙网络，这可解释为沉淀溶解反应作用的结果. 在120°的网络中，Pb2+生成的沉淀量相对较少，从而促成了最终浓度的增高，形成一种拮抗作用. 一般实际工程应用中的目标都是尽量降低地下水中的污染物浓度，对于保守性溶质而言，当地下水总体流向和整体渗透场最大渗透主轴方向不平行时，污染物迁移的范围就能够被压缩，但在某些特定的反应条件下，平行于渗透主轴方向的地下水流向反而能促成更为充分的沉淀反应，由此可降低污染的浓度. 对于Pb2+而言，不同方向的裂隙网络在最终下游断面上展示了120°走向的裂隙网络超过90°网络的平均浓度. 另一方面，对比这些网络的下游断面通量，如图 6d所示，90°、60°、120°对应裂隙网络在最终1 000 d左右的最终通量非常接近，这是因为它们在最终断面上的流速分布差异小，而75°的裂隙网络由于污染源没有连通主要裂隙，且大气降水在持续稀释Pb2+浓度，导致其最终通量趋于零；60°和120°裂隙网络在Pb2+浓度上的对称性虽不明显，但在最终的通量上两者没有明显的差异.
3.3 裂隙立方定律修正对运移结果的影响
由公式（20）可知，对立方定律是否进行修正以及如何修正，其实质就是如何刻画同一裂隙网络的渗透性，由此会进一步影响溶质的运移结果. 如表 4所示，随着修正系数的增大，裂隙网络的渗透率增强，其中的水流速度加快，由对流迁移的各组分质量随之增多，造成下游断面的溶质浓度降低（图 7a），而其对应的总通量增大（图 7b）. 如图 7所示，在高渗透率网络情景下，最终下游断面溶质的平均浓度更低，但整个断面的通量则更高. 同时，图 7c显示反应性组分（如Pb2+）与保守性组分（如Cl-）具有类似的变化规律.
图 7 不同渗透率裂隙网络情景下不同组分在下游断面的穿透曲线a. Cl-浓度；b. Cl-通量；c. Pb2+浓度Fig. 7. Breakthrough curve of different components in the downstream section under different permeability fracture network scenarios若比较同一组分的溶质等值线锋面，可发现渗透率越低，同一组分的浓度等值锋面距离初始污染源越远. 如图 8所示，通过比较第1 000 d切面上的Cl-浓度分布，可知渗透率同时亦影响了岩石骨架中的溶质弥散：渗透率越低的网络，其骨架中的弥散越明显，进而削减其最终在下游断面的溶质通量.
图 8 不同渗透率裂隙网络情景下同一切面位置上对应的Cl-浓度分布Fig. 8. Distribution of Cl- concentration at the same slice location under different permeability fracture network scenarios对立方定律修正实质就是考虑粗糙度对渗透率的影响. 对于通量而言，最终流速的影响往往强于浓度的影响，而从数值上来说，流速在裂隙和岩石骨架中的变化会横跨若干数量级，比浓度的变化大得多. 实际场地中裂隙并不光滑，理应考虑裂隙的粗糙对立方定律进行修正. 至于如何修正本文不予讨论，有待于后续结合场地监测数据进一步研究，通过拟合得到合适的修正系数来确刻画裂隙网络中的流场和浓度场，以准确评估污染物的总通量.
4. 结论
（1）通过构建不同情景的DFN并映射至EPM，使用PFLOTRAN代码模拟了一定源强污染物连续释放情况下各组分在裂隙网络中的迁移规律，并量化评估其对应特定断面的组分通量. 从均质多孔介质到基于DFN映射至EPM，由二维剖面到三维裂隙网络，本文为真实裂隙岩体场地的复杂条件下多组分反应性污染物运移过程的数值模拟提供一种视角.
（2）算例结果表明，裂隙发育的密度和方向均影响着裂隙网络的连通性，进而影响裂隙水流中不同组分的迁移规律和通量评估，其中保守性组分与不同反应性组分的迁移行为有明显差异. 裂隙介质的岩石骨架对污染物迁移有一定的吸收和阻滞作用，且在pH、碳酸盐岩及硫酸根离子的作用下，某些重金属污染物（如文中Pb2+）最终会在特定位置被沉淀固定下来.
（3）三维裂隙网络的特定几何拓扑结构对模拟污染运移行为和通量评估具有重要影响. 本文方法并没有对裂隙及骨架进行精细网格剖分以规避计算负荷. 区域尺度上结合实际地质背景对裂隙网络做出更加合理的刻画（Bense et al.，2013）有待于后续进一步探索如何有效处理非均质性和各项异性带来的复杂度是能够合理使用模型进行评估的重要前提.
-
图 4 高密度裂隙网络中不同组分运移结果及其空间分布
a.对流速度场分布及保守性Cl-浓度空间分布；b. Cl-、Fe2+浓度以及菱铁矿饱和指数空间分布；c.Pb2+浓度及碳酸铅、硫酸铅、软石膏饱和指数的空间分布
Fig. 4. Solute transport results and their spatial distributions of different components in the high-density fracture network
图 5 不同密度裂隙网络情景下不同组分在下游断面的浓度穿透曲线
a. Cl-平均浓度；b. Fe2+最大浓度
Fig. 5. Breakthrough curves of different components in the downstream section under different density fracture network scenarios
图 6 不同走向裂隙网络情景下不同组分在下游断面的穿透曲线
a. Cl-平均浓度；b. Fe2+平均浓度；c. Pb2+平均浓度；d. Pb2+通量
Fig. 6. Breakthrough curves of different components in the downstream section under different trends fracture network scenarios
图 7 不同渗透率裂隙网络情景下不同组分在下游断面的穿透曲线
a. Cl-浓度；b. Cl-通量；c. Pb2+浓度
Fig. 7. Breakthrough curve of different components in the downstream section under different permeability fracture network scenarios
图 8 不同渗透率裂隙网络情景下同一切面位置上对应的Cl-浓度分布
Fig. 8. Distribution of Cl- concentration at the same slice location under different permeability fracture network scenarios
表 1 溶质运移模型初始及边界各组分浓度（单位：mol$ \cdot $L-1, pH除外）
Table 1. Initial and boundary condition of the solute transport model (unit: mol$ \cdot $L-1, except pH)
主要离子组分 背景值 源强组分 降水组分 H+ (pH) 9.2 4.0 7.0 Al3+ 1.00×10-20 4.30×10-3 1.30×10-7 Fe3+ 1.00×10-22 2.00×10-7 2.30×10-8 Fe2+ 4.20×10-5 3.06×10-2 5.40×10-5 Mn2+ 2.90×10-5 7.84×10-3 4.70×10-5 SO42- 5.30×10-3 5.00×10-2 7.50×10-3 CO32- 2.10×10-5 4.92×10-4 3.90×10-3 Cl- 1.00×10-3 1.58×10-2 1.00×10-3 K+ 6.60×10-5 8.14×10-4 6.70×10-5 Na+ 1.30×10-3 1.39×10-3 1.30×10-3 Mg2+ 1.50×10-3 9.69×10-4 1.90×10-3 Ca2+ 5.20×10-3 1.08×10-2 6.09×10-3 Pb2+ 1.00×10-8 1.52×10-5 1.00×10-20 
下载: 导出CSV
表 2 不同密度连通裂隙网络的参数
Table 2. Parameters of connected fracture networks with different density
裂隙密度 目标裂隙数N/(条) 连通裂隙数n/(条) P32 低密度 250 54 0.16 中密度 500 277 0.70 高密度 1 000 511 1.12
下载: 导出CSV
表 3 中密度（N=500）条件下生成不同走向的裂隙网络特征
Table 3. Fracture network features with different trends generated under the condition of medium density (N=500)
不同情景 走向φ(°) 连通裂隙数n(条) P32 Case 0（基准方向） 90 277 0.70 Case 1 60 192 0.51 Case 2 75 92 0.24 Case 3 105 114 0.32 Case 4 120 256 0.66
下载: 导出CSV
表 4 不同渗透率裂隙网络设定
Table 4. Fracture networks settings with different permeability
渗透特性 修正系数 渗透率取值范围(m2) $ \overline{{\boldsymbol{k}}_{\boldsymbol{x}\boldsymbol{x}}} $ $ \overline{{\boldsymbol{k}}_{\boldsymbol{y}\boldsymbol{y}}} $ $ \overline{{\boldsymbol{k}}_{\boldsymbol{z}\boldsymbol{z}}} $ 低渗透率 0.24 [2.95×10-13, 3.58×10-12] [5.32×10-17, 1.54×10-12] [3.18×10-13, 3.47×10-12] 中渗透率 0.48 [5.78×10-13, 7.02×10-12] [1.04×10-16, 3.02×10-12] [6.22×10-13, 6.79×10-12] 高渗透率 1.00 [1.21×10-12, 1.47×10-11] [2.19×10-16, 6.35×10-12] [1.31×10-12, 1.43×10-11]
下载: 导出CSV
表 5 化学反应及其平衡常数
Table 5. Chemical reactions and their equilibrium constants
反应表达式 平衡常数lgK $ \mathrm{O}{\mathrm{H}}^{-}+{\mathrm{H}}^{+}\leftrightarrow {\mathrm{H}}_{2}\mathrm{O} $ 13.995 1 $ {\mathrm{H}}^{+}+\mathrm{C}{\mathrm{O}}_{3}^{2-}\leftrightarrow \mathrm{H}\mathrm{C}{\mathrm{O}}_{3}^{-} $ 10.328 8 $ \mathrm{C}{\mathrm{O}}_{2}\left(\mathrm{a}\mathrm{q}\right)+{\mathrm{H}}_{2}\mathrm{O}\leftrightarrow {\mathrm{H}}^{+}+\mathrm{H}\mathrm{C}{\mathrm{O}}_{3}^{-} $ -6.344 7 $ \mathrm{C}\mathrm{a}\mathrm{C}{\mathrm{O}}_{3}+{\mathrm{H}}^{+}\leftrightarrow \mathrm{C}{\mathrm{a}}^{2+}+\mathrm{H}\mathrm{C}{\mathrm{O}}_{3}^{-} $ 1.848 7 $ \mathrm{F}\mathrm{e}\mathrm{C}{\mathrm{O}}_{3}+{\mathrm{H}}^{+}\leftrightarrow \mathrm{F}{\mathrm{e}}^{2+}+\mathrm{H}\mathrm{C}{\mathrm{O}}_{3}^{-} $ -0.192 0 $ \mathrm{A}\mathrm{l}{\left(\mathrm{O}\mathrm{H}\right)}_{3}+3{\mathrm{H}}^{+}\leftrightarrow \mathrm{A}{\mathrm{l}}^{3+}+3{\mathrm{H}}_{2}\mathrm{O} $ 7.756 0 $ \mathrm{F}\mathrm{e}{\left(\mathrm{O}\mathrm{H}\right)}_{3}+3{\mathrm{H}}^{+}\leftrightarrow \mathrm{F}{\mathrm{e}}^{3+}+3{\mathrm{H}}_{2}\mathrm{O} $ 4.896 0 $ \mathrm{C}\mathrm{a}\mathrm{S}{\mathrm{O}}_{4}\cdot 2{\mathrm{H}}_{2}\mathrm{O}\leftrightarrow \mathrm{C}{\mathrm{a}}^{2+}+\mathrm{S}{\mathrm{O}}_{4}^{2-}+2{\mathrm{H}}_{2}\mathrm{O} $ -4.482 3 $ \mathrm{P}\mathrm{b}\mathrm{S}{\mathrm{O}}_{4}\leftrightarrow \mathrm{P}{\mathrm{b}}^{2+}+\mathrm{S}{\mathrm{O}}_{4}^{2-} $ -7.852 7 $ \mathrm{P}\mathrm{b}\mathrm{C}{\mathrm{O}}_{3}+{\mathrm{H}}^{+}\leftrightarrow \mathrm{P}{\mathrm{b}}^{2+}+\mathrm{H}\mathrm{C}{\mathrm{O}}_{3}^{-} $ -3.209 1
下载: 导出CSV
-
Barton, N., Bandis, S., Bakhtar, K., 1985. Strength, Deformation and Conductivity Coupling of Rock Joints. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 22(3): 121-140. https://doi.org/10.1016/0148-9062(85)93227-9 Bense, V. F., Gleeson, T., Loveless, S. E., et al., 2013. Fault Zone Hydrogeology. Earth-Science Reviews, 127: 171-192. https://doi.org/10.1016/j.earscirev.2013.09.008 Chen, H. F., Zhang, F. W., He, Y., et al., 2016. Geological and Geomorphologic Settings Acting as the Controlling Factors and Indicators for Karst Systems. Hydrogeology & Engineering Geology, 43(5): 42-47(in Chinese with English abstract). de Dreuzy, J. R., Méheust, Y., Pichot, G., 2012. Influence of Fracture Scale Heterogeneity on the Flow Properties of Three-Dimensional Discrete Fracture Networks (DFN). Journal of Geophysical Research: Solid Earth, 117(B11) https://doi.org/10.1029/2012JB009461 Dershowitz, W. S., Einstein, H. H., 1988. Characterizing Rock Joint Geometry with Joint System Models. Rock Mechanics and Rock Engineering, 21(1): 21-51. https://doi.org/10.1007/BF01019674 Gan, L., Ma, H. Y., Shen, Z. Z., 2021. Roughness Characterization of Concave Fracture Surface and Coefficient Fitting of Modified Cubic Law. Journal of Hydraulic Engineering, 52(4): 420-431(in Chinese with English abstract). Hammond, G. E., Lichtner, P. C., 2010. Field-Scale Model for the Natural Attenuation of Uranium at the Hanford 300 Area Using High-Performance Computing. Water Resources Research, 46(9). https://doi.org/10.1029/2009WR008819 Hu, Y., Xu, W., Zhan, L., et al., 2022. Modeling of Solute Transport in a Fracture-Matrix System with a Three-Dimensional Discrete Fracture Network. Journal of Hydrology, 605: 127333. https://doi.org/10.1016/j.jhydrol.2021.127333 Hyman, J. D., Aldrich, G., Viswanathan, H., et al., 2016. Fracture Size and Transmissivity Correlations: Implications for Transport Simulations in Sparse Three-Dimensional Discrete Fracture Networks Following a Truncated Power Law Distribution of Fracture Size. Water Resources Research, 52(8): 6472-6489. https://doi.org/10.1002/2016WR018806 Hyman, J. D., Gable, C. W., Painter, S. L., et al., 2014. Conforming Delaunay Triangulation of Stochastically Generated Three Dimensional Discrete Fracture Networks: A Feature Rejection Algorithm for Meshing Strategy. SIAM Journal on Scientific Computing, 36(4): A1871-A1894. https://doi.org/10.1137/130942541 Hyman, J. D., Karra, S., Makedonska, N., et al., 2015. dfnWorks: A Discrete Fracture Network Framework for Modeling Subsurface Flow and Transport. Computers & Geosciences, 84: 10-19. https://doi.org/10.1016/j.cageo.2015.08.001 Jordan, G., Rammensee, W., 1998. Dissolution Rates of Calcite (1014) Obtained by Scanning Force Microscopy: Microtopography-Based Dissolution Kinetics on Surfaces with Anisotropic Step Velocities. Geochimica et Cosmochimica Acta, 62(6): 941-947. https://doi.org/10.1016/S0016-7037(98)00030-1 Lei, Q., Latham, J. P., Tsang, C. F., 2017. The Use of Discrete Fracture Networks for Modelling Coupled Geomechanical and Hydrological Behaviour of Fractured Rocks. Computers and Geotechnics, 85: 151-176. https://doi.org/10.1016/j.compgeo.2016.12.024 Li, X., Li, D., Xu, Y., et al., 2020. A DFN based 3D numerical Approach for Modeling Coupled Groundwater Flow and Solute Transport in Fractured Rock Mass. International Journal of Heat and Mass Transfer, 149: 119179. https://doi.org/10.1016/j.ijheatmasstransfer.2019.119179 Li, Y. M., Wen, Z., 2020. Impacts of Non-Darcian Flow in the Fracture on Flow Field and Solute Plumes in a Fracture-Aquifer System. Earth Science, 45(2): 693-700(in Chinese with English abstract). Liang, M., 2010. Numerical Simulation Study of Flow in Single Rough Fractures Based on Fluent(Dissertation). Hefei University of Technology, Hefei(in Chinese with English abstract). Lichtner, P. C., 2016. Kinetic Rate Laws Invariant to Scaling the Mineral Formula Unit. American Journal of Science, 316(5): 437. https://doi.org/10.2475/05.2016.02 Long, J. C. S., Karasaki, K., Davey, A., et al., 1990. Preliminary Prediction of Inflow into the D-Holes at the Stripa Mine. Lawrence Berkeley National Laboratory. LBNL Report #: LBL-27182. https://escholarship.org/uc/item/7mb9q37h Makedonska, N., Hyman, J. D., Karra, S., et al., 2016. Evaluating the Effect of Internal Aperture Variability on Transport in Kilometer Scale Discrete Fracture Networks. Advances in Water Resources, 94: 486-497. https://doi.org/10.1016/j.advwatres.2016.06.010 Makedonska, N., Painter, S. L., Bui, Q. M., et al., 2015. Particle Tracking Approach for Transport in Three-Dimensional Discrete Fracture Networks. Computational Geosciences, 19(5): 1123-1137. https://doi.org/10.1007/s10596-015-9525-4 Molson, J., Aubertin, M., Bussière, B., et al., 2008. Geochemical Transport Modelling of Drainage from Experimental Mine Tailings Cells Covered by Capillary Barriers. Applied Geochemistry, 23(1): 1-24. https://doi.org/10.1016/j.apgeochem.2007.08.004 Molson, J., Aubertin, M., Bussière, B., 2012. Reactive Transport Modelling of Acid Mine Drainage within Discretely Fractured Porous Media: Plume Evolution from a Surface Source Zone. Environmental Modelling & Software, 38: 259-270. https://doi.org/10.1016/j.envsoft.2012.06.010 Pandey, S., Rajaram, H., 2016. Modeling the Influence of Preferential Flow on the Spatial Variability and Time-Dependence of Mineral Weathering Rates. Water Resources Research, 52(12): 9344-9366. https://doi.org/10.1002/2016WR019026 Peters, R. R., Klavetter, E. A., 1988. A Continuum Model for Water Movement in an Unsaturated Fractured Rock Mass. Water Resources Research, 24(3): 416-430. https://doi.org/10.1029/WR024i003p00416 Romano, V., Bigi, S., Carnevale, F., et al., 2020. Hydraulic Characterization of a Fault Zone from Fracture Distribution. Journal of Structural Geology, 135: 104036. https://doi.org/10.1016/j.jsg.2020.104036 Sweeney, M. R., Gable, C. W., Karra, S., et al., 2020. Upscaled Discrete Fracture Matrix Model (UDFM): an Octree-Refined Continuum Representation of Fractured Porous Media. Computational Geosciences, 24(1): 293-310. https://doi.org/10.1007/s10596-019-09921-9 Trinchero, P., Iraola, A., Bruines, P., et al., 2021. Water-Mineral Reactions in a Translated Single Realistic Fracture: Consequences for Contaminant Uptake by Matrix Diffusion. Water Resources Research, 57(10): e2021WR030442. https://doi.org/10.1029/2021WR030442 Tsang, Y. W., Tsang, C. F., 1987. Channel Model of Flow through Fractured Media. Water Resources Research, 23(3): 467-479. https://doi.org/10.1029/WR023i003p00467 Walter, A. L., Frind, E. O., Blowes, D. W., et al., 1994. Modeling of Multicomponent Reactive Transport in Groundwater: 2. Metal Mobility in Aquifers Impacted by Acidic Mine Tailings Discharge. Water Resources Research, 30(11): 3149-3158. https://doi.org/10.1029/94WR00954 Wang, W., Fu, H., Xing, L. X., et al., 2021. Crack Propagation Behavior of Carbonatite Geothermal Reservoir Rock Mass Based on Extended Finite Element Method. Earth Science, 46(10): 3509-3519(in Chinese with English abstract). White, A. F., Brantley, S. L., 2003. The Effect of Time on the Weathering of Silicate Minerals: Why do Weathering Rates Differ in the Laboratory and Field? Chemical Geology, 202(3): 479-506. https://doi.org/10.1016/j.chemgeo.2003.03.001 Wilson, C. R., Witherspoon, P. A., Long, J. C. S., et al., 1983. Large-Scale Hydraulic Conductivity Measurements in Fractured Granite. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 20(6): 269-276. https://doi.org/10.1016/0148-9062(83)90596-X Yang, B., Wang, H., Wang, B., et al., 2021. Digital Quantification of Fracture in Full-Scale Rock Using Micro-CT Images: A Fracturing Experiment with N2 and CO2. Journal of Petroleum Science and Engineering, 196: 107682. https://doi.org/10.1016/j.petrol.2020.107682 Zhong, Q. M., Chen, J. S., Chen, L., 2006. Symmetry Certification of Permea Bility Tensor and Derivation of Principal Permea Bility of Fractured Rock Mass. Chinese Journal of Rock Mechanics And Engineering, (S1): 2997-3002(in Chinese with English abstract). Zimmerman, R. W., Yeo, I. W., 2000. Fluid Flow in Rock Fractures: From the Navier-Stokes Equations to the Cubic Law. Dynamics of Fluids in Fractured Rock: 213-224. https://doi.org/10.1029/GM122p0213 陈宏峰, 张发旺, 何愿, 等, 2016. 地质与地貌条件对岩溶系统的控制与指示. 水文地质工程地质, 43(5): 42-47. https://www.cnki.com.cn/Article/CJFDTOTAL-SWDG201605006.htm 甘磊, 马洪影, 沈振中, 2021. 下凹形态裂隙面粗糙程度表征及立方定律修正系数拟合. 水利学报, 52(4): 420-431. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB202104005.htm 梁敏, 2010. 基于Fluent的粗糙单裂隙水流数值模拟研究(硕士学位论文). 合肥: 合肥工业大学. 李一鸣, 文章, 2020. 非达西裂隙流对渗透性基岩中流场及溶质羽的影响. 地球科学, 45(2): 693-700. doi: 10.3799/dqkx.2018.345 王伟, 付豪, 邢林啸, 2021. 基于扩展有限元法的碳酸盐岩地热储层岩体裂缝扩展行为. 地球科学, 46(10): 3509-3519. doi: 10.3799/dqkx.2021.005 钟启明, 陈建生, 陈亮, 2006. 裂隙岩体渗透张量的对称性证明及主渗透性推导. 岩石力学与工程学报, (S1): 2997-3002. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2006S1060.htm 期刊类型引用(1)
1. 刘殿广，杨蕴，袁奕梁，窦智，吴剑锋，吴吉春. 基于升尺度离散裂隙—基质模型的裂隙水反应运移数值模拟与软件开发. 高校地质学报. 2025(01): 14-23 . 
百度学术其他类型引用(1)
-
下载:
百度学术
点击查看大图
计量
- 文章访问数: 490
- HTML全文浏览量: 203
- PDF下载量: 78
- 被引次数: 2
