一种用于表层岩溶区高分辨率裂缝网络的情感约束简化方法

《Journal of Hydrology: Regional Studies》：A flow constrained simplification method for high-resolution fracture networks in epikarst zone

　　高风军|陈曦|张志才|梁盈春|刘威尔汉
河海大学水文学与水资源学院，南京210098，中国

**摘要**
**研究区域**
本研究在贵州省普定县进行。

**研究重点**
喀斯特地区广泛存在裂隙含水层。这些裂隙网络通常具有复杂的结构，这使得流动和溶质传输的模拟在计算上非常耗费资源且数值上不稳定。这严重限制了我们对喀斯特裂隙系统中不确定性的量化和表征能力。本研究提出了一种简化方法，用于从现场露头中提取的高分辨率裂隙网络。通过使用结构复杂性和尺度不同的天然裂隙网络来评估所提出方法的性能和适用性。

**该地区的新水文见解**
对四个天然裂隙网络的数值模拟表明，尽管高分辨率裂隙模型提供了详细的几何信息，但并非所有几何细节都对宏观水力响应具有控制性影响。过度的几何细化在提高水力表示方面效果有限，同时显著增加了网格密度和自由度，从而降低了计算效率。所提出的流动约束简化方法在有效保留主要流动路径和等效渗透性的同时，减少了局部结构复杂性。溶质传输模拟进一步证明，在适当选择的阈值下，突破曲线偏差保持在可接受的范围内，并且明显小于使用无约束几何简化方法得到的偏差，尤其是对于结构复杂的网络。所提出的方法提供了一个基于物理原理且计算效率高的框架，可以在不牺牲关键水文过程的情况下降低模型复杂性。

**1. 引言**
喀斯特地貌占地球陆地表面的大约20%（Bakalowicz, 2005）。作为喀斯特表面附近的薄层或“皮肤”，表层喀斯特具有丰富的裂隙网络特征，这些特征是由构造应力、溶解作用和生物活动形成的。裂隙网络为地下水流动和溶质传输提供了主要通道（Berkowitz, 1995; Gan和Elsworth, 2016; Sanderson和Nixon, 2018; Tokunaga和Wan, 2001; Wang和Cardenas, 2014）。然而，裂隙结构和网络拓扑的固有复杂性，加上其中流动过程的非线性特性，使得准确高效地进行地下水流动和溶质传输的数值建模成为科学挑战。

**离散裂隙矩阵（DFM）模型**
DFM模型是一种广泛用于高保真模拟裂隙多孔介质中流动和传输的方法，其中每个裂隙都被明确表示并赋予特定的几何和水力属性（Berre等, 2019; Lei等, 2017; Ngo等, 2017）。在DFM中，裂隙通常被表示为网格中的低维对象（例如，在二维中为线条，在三维中为圆盘），这样可以捕捉裂隙与周围基质之间的交叉流动平衡（Hyman, 2020; Zhang等, 2022）。DFM高度依赖于裂隙网络的几何属性。近年来，地球物探和图像处理技术的进步使得裂隙网络的高分辨率表征变得越来越可行（Bemis等, 2014; Bisdom等, 2017; Chen等, 2021; Pan等, 2019; Viswanathan等, 2022）。与统计生成的模型或低分辨率的大规模模型（从米到公里）不同，这些技术可以轻松地从露头剖面中捕获高分辨率的裂隙网络，同时保留详细的几何属性，如开口度、方向和间距，以及裂隙的完整拓扑连接性。这些露头特征为模拟现实裂隙系统中的流动和溶质传输提供了宝贵的数据集。

DFM模型的高计算成本源于共形网格的使用和裂隙网络的内在复杂性。成千上万个明确表示的裂隙需要精细的空间离散化，特别是在低角度交点、密集排列的裂隙或高度曲折的路径处。这些几何特征需要局部细化的网格以确保数值稳定性和准确性，这显著增加了计算负担。本质上，捕捉详细的裂隙几何形状不可避免地导致自由度的大幅增加（De Hoop等, 2022; Zou等, 2017）。此外，裂隙网络属性的不确定性量化通常需要数千次模拟来捕捉与单个裂隙的位置、长度、方向和水力属性相关的变异性（Lei等, 2020; Sweeney等, 2023; Xue等, 2023）。因此，有必要在保持DFM准确性的同时减少计算时间和成本。

为了减轻裂隙网络建模的高计算成本，提出了各种简化策略，以在保留关键水力特征的同时降低裂隙网络复杂性并提高计算效率。几何简化通过局部调整裂隙属性（如长度、节点位置或方向）来改善生成网格的质量（Anquez等, 2019; Karimi-Fard和Durlofsky, 2016; Mustapha和Dimitrakopoulos, 2011; Mustapha和Mustapha, 2007），但它对模型规模的减小效果有限，因为裂隙的总数通常会保留。相比之下，拓扑简化通过评估裂隙的相对水力重要性来保留主要流动路径，同时移除或合并贡献较小的裂隙，从而在不破坏主要流动结构的情况下实现显著的模型简化。这些方法包括路径搜索算法（例如最短路径）、基于拓扑的方法（例如中心性或节点度），以及用于识别骨干子网络的数据驱动技术（Valera等, 2018; Viswanathan等, 2018; Zhu等, 2021）。虽然这些简化方法可以显著减小模型规模并提高计算效率，但它们通常依赖于几何或网络拓扑指标，而不是直接考虑流场。由于在指定边界条件下典型的DFM中的流动行为是唯一且确定性的，任何对网络拓扑的修改都可能影响这种独特的流场（Wellman等, 2009）。因此，裂隙的移除或合并可能导致流动和溶质传输的预测偏差（Viswanathan等, 2018; Zhang等, 2024）。正如Zhu等（2021）指出的，如果裂隙网络具有分形空间密度分布并以小裂隙为主，那么基于最小阻力路径的简化方法是不足的。然而，这些方法主要依赖于几何和拓扑标准来选择要移除的裂隙，而没有明确评估哪些特定的删除会改变流动路径，哪些会保持原始网络的全局水力行为和关键流动结构。特别是，如何从露头剖面图像中简化真实的裂隙网络以降低计算成本，同时保留关键的水力行为，这一问题尚未得到很好的解决。

本研究旨在开发一种适用于喀斯特地区复杂裂隙网络的简化方法，特别是适用于从露头中提取的高分辨率网络。该方法可以在可接受的精度范围内显著提高计算效率，从而为水文地质参数不确定性的定量分析和大规模情景模拟提供可行的数值工具。我们的方法将裂隙网络的拓扑结构与流场信息相结合。通过流动约束，它识别并保留控制宏观流动的关键裂隙骨架，同时系统地补偿由几何调整引起的水力变化。

本文的其余部分组织如下。第2节描述了从现场露头中提取高分辨率裂隙网络的过程，介绍了所提出的流动约束简化算法的理论框架，并提出了流动和溶质传输建模的控制方程。第3节评估了简化网络的结构演变，并从等效渗透性、突破行为和计算效率等方面定量评估了它们的水力和传输性能。最后，第4节将该方法应用于站点规模的裂隙网络，并详细讨论了其适用性、局限性和鲁棒性。

**2. 材料与方法**
2.1. 选定剖面露头中裂隙网络的特征
2.1.1. 从剖面露头图像中提取裂隙网络
在中国贵州省普定县的喀斯特地区选择了三个典型的剖面露头（图1）。由于该地区位于亚热带湿润气候区，物理和化学风化作用强烈，这三个剖面都具有丰富的裂隙。在地质形成方面，所有剖面中的主要岩性是石灰岩，属于中三叠世的观音群（T2g），包括图1a和1b中的T2g2–2亚组和图1c中的T2g2–1。

从每个剖面中选择1×1米的露头部分，并使用DSLR相机（SONY α6000）拍摄剖面图像。然后，使用LIDAR（光检测和测距）扫描数据对图像进行几何校正，以确保准确提取裂隙的几何特征（Lee等, 2022; Marques等, 2020）。捕获了多种分辨率的图像，以增强对更细小裂隙的检测，像素尺寸从0.1毫米到1.7毫米不等。手动解释露头图像中的裂隙以生成二值图像，然后用于提取多段线和随后的裂隙网络简化。

采用了一种基于像素的裂隙提取方法，包括以下步骤（图2）：首先，使用MATLAB中的‘bwskel’函数对二值裂隙图像进行骨架化，得到单像素宽度的中心线（图2b）。随后使用基于链码的骨干提取方法跟踪裂隙轨迹（图2c），该方法可以检测八个方向上的像素连通性并生成连续的裂隙坐标（Chen等, 2021）。每个裂隙的开口度定义为从中心线到两侧边缘的最短距离之和的平均值。从现场测量得到的裂隙开口度与从图像中得到的开口度之间的相对误差不超过5%。露头图像的手动解释可能会引入小的位置误差，因此，视觉上似乎相交的裂隙可能在端点上不完全重合（Bisdom等, 2017; Zhu等, 2022）。为了纠正这些不一致性并保持几何连通性，检查了所有裂隙端点，并将距离在3像素以内的端点合并。

**2.2. 裂隙网络的拓扑和几何**
裂隙网络可以通过拓扑属性（Healy等, 2017; Sanderson和Nixon, 2018）来定量描述，例如顶点、边以及网络图中顶点的度。在本研究中，顶点定义为裂隙的交点或端点，而边表示连接两个顶点的裂隙段。顶点的度对应于连接到该顶点的边的数量。统计上，这些拓扑参数可以通过度分布的均值（μ(b)）和方差（σ(b)）以及香农熵（H）来表示，它们表明了通过图/网络的流动和传输能力。注意：H=?∑p(c)log2pc?1，其中pc是图中的度中心性分布（Hyman, 2020）；裂隙密度（P）定义为单位面积内的总裂隙长度。

从图1中的三个剖面露头中提取的裂隙网络结构显示在图3中，它们的密度和拓扑参数总结在表1中。网络A（图3a）的裂隙密度最低（10.07米/平方米），顶点数量最少（145个）。其结构相对简单，具有明确的优先流动路径。平均节点度为2.97，方差较大，而香农熵最小（5.55），表明拓扑简单，水力连通性集中在少数主导通道中。

网络B（图3b）的裂隙密度较大，为25.72米/平方米，有656个顶点。几何上，裂隙较短且较窄，大约90%的裂隙开口度在0.5–2毫米之间。网络由更多的水平和垂直裂隙交点组成，而垂直裂隙的优先路径较少。平均节点度为3.42，方差相对较小，表明连通性更加均匀。香农熵达到7.56，显著高于网络A，反映了网络结构的复杂性和流动路径的多样性。

网络C（图3c）的裂隙密度最高（28.09米/平方米），顶点数量最多（1021个）。其结构特征是密集的短裂缝群集，形成了一个高度复杂的网络。在入口和出口之间存在多条平行或替代的流动路径，导致流体分布更加分散。平均节点度为3.34，方差较大，表明节点连接性存在显著异质性。香农熵最高（8.20），表明拓扑结构最为复杂，流动路径的多样性也最大。2.2. 裂缝网络简化2.2.1. 裂缝网络简化算法从露头图像中提取的数字化裂缝网络通常包含高密度的节点、众多的短裂缝段和小角度交叉点，这些因素大大增加了网络的拓扑复杂性并降低了计算效率。为了降低原始裂缝网络模型的复杂性，本研究提出了一种基于流动约束的裂缝网络简化方法。该方法综合考虑了网络拓扑结构和稳态水力条件，以提取主要的裂缝骨架，同时保留宏观水力响应。在提出的裂缝网络简化方法中，假设石灰岩基质的渗透率远低于裂缝的渗透率，并且地下水流动主要由裂缝控制（刘等人，2025；刘等人，2024）。因此，可以忽略基质流动的贡献。在此框架内，流动约束的核心在于保持局部水力梯度的稳定性。任何简化过程中的拓扑修改（例如，裂缝去除、合并或端点调整）都不得改变相邻段的流动方向。如果发生流动反向，则表明局部水力梯度受到干扰，可能会影响连接到该修改段的所有路径的流动；因此，这种操作是被禁止的。当没有发生流动反向时，可以通过调整裂缝开口来补偿由于几何变化引起的水量不平衡。在实际操作中，所有简化操作都严格受到原始网络稳态流场的约束。每次拓扑修改仅限于局部区域，而未修改段的流动方向和大小保持不变，从而保护关键流动路径，并确保简化后的网络整体水力行为与原始网络一致。如图4所示，简化过程包括两个主要步骤。第一步是通过对裂缝段进行几何平滑处理来减少节点密度和裂缝迂曲度。第二步侧重于拓扑调整，包括删除短裂缝段和消除锐角，目的是在保持网络连通性的同时降低结构复杂性。下载：下载高分辨率图像（374KB）下载：下载全尺寸图像图4. 从原始裂缝网络到简化裂缝网络的简化工作流程。lv表示节点间距阈值：它控制裂缝段的平滑处理，并用作删除短段的长度阈值。φ表示消除锐角的角度阈值。裂缝网络简化步骤的详细信息如下：步骤1：几何平滑：使用预设的节点间距阈值lv对每个裂缝段进行几何平滑处理，以减少迂曲度和节点密度（图5b）。长度小于lv的段被简化为直线，而长度大于lv的段则使用三次样条曲线进行平滑处理。平滑后，调整每个段的开口大小以补偿由于段长度缩短引起的水力变化。随着lv的增加，节点数量显著减少，几何形状变得更加平滑。需要强调的是，这一步仅修改裂缝段的几何形状，而不改变它们的连通性；因此，网络拓扑保持不变。下载：下载高分辨率图像（341KB）下载：下载全尺寸图像图5. 从原始网络到简化裂缝网络的简化步骤，其中实心框和虚线框分别代表处理前后的状态。步骤2：基于流动约束的拓扑调整。(1)删除短段：为了最小化拓扑扰动并保持原始网络结构，选择短于阈值lv的最短裂缝段作为删除候选，其端点分别用A和B表示（图5d）。首先，将原本连接到端点A的裂缝段重新连接到端点B，并使用原始网络的水力头分布来评估受影响段的流动方向。如果在任何修改后的段中检测到流动反向，则认为这种调整是不可接受的。然后尝试另一种调整方法，即将原本连接到端点B的段重新连接到端点A，并进行相同的流动一致性检查。如果两种调整都导致流动反向，则拒绝删除目标裂缝段。算法随后继续处理下一条最短的裂缝段，并重复上述调整和一致性检查。否则，当流动方向保持一致时，调整段的开口大小以补偿由于结构修改引起的局部流动变化，确保流动平衡。(2)消除锐角：删除裂缝段后，必须检查网络中是否存在低于角度阈值φ的锐角，无论是预先存在的还是由于结构调整新产生的，因为这些角度会降低网格质量（图5f）。遍历所有节点，并计算每个节点处相邻边之间的角度。如果角度小于阈值φ，则合并属于不同裂缝段的边以消除锐角，同时通过移除节点来伸直属于同一裂缝段的边，以消除过度弯曲。所有锐角调整都需经过流动方向一致性和开口大小校正，以保持质量守恒。步骤2的子步骤（1）和（2）重复进行，直到所有长度小于lv的裂缝段要么被删除，要么由于删除后会导致相邻段流动反向而被认定为不可移除，此时迭代过程结束。为了保持网络的整体连通性，在简化过程中保留所有直接连接到模型边界的裂缝段。每次拓扑调整后，更新受影响段的水力信息（流量和水头），以避免重复进行全局流动计算，从而降低计算成本。锐角阈值的选取是基于网格数量和网格质量之间的权衡（补充图1）。当角度低于20°时，网格数量迅速增加且网格质量下降，而对于较大的角度，变化较小。为了在减少拓扑扰动的同时提高计算效率，采用了20°的交叉角阈值。虽然较小的阈值更有利于保留原始网络的详细几何信息，但它们伴随的网格密度增加和网格质量下降会大幅降低计算效率。2.2.2. 裂缝开口校正简化过程改变了长度和连通性等几何参数。为了保持等效的流动条件，通过调整受影响裂缝的开口大小来校正水力传导性。如图6a和6b所示，短段删除和/或端点位置的变化会改变流动长度，从而改变流量，这可以通过调整裂缝开口来补偿。此外，在删除裂缝段后，端点合并或几何调整可能导致相邻裂缝部分或完全重叠，这些裂缝合并成一个等效裂缝，并应用开口校正，如图6c所示。下载：下载高分辨率图像（95KB）下载：下载全尺寸图像图6. 网络简化中的开口校正：(a) 节点删除；(b) 删除裂缝段后的端点位移，虚线表示被删除的裂缝；(c) 重叠的裂缝段。线条颜色区分不同的裂缝段。简化后，修改后的裂缝段应保持与原始条件相同的渗透能力。根据任何段中的流量守恒原理，三种网络简化方法的校正可以表示为：(1)∑i=1nQori,i=Kdsimh1?h2Lsim其中Qori,i是通过第i个原始裂缝段的流量（包含方向信息），K是水力传导率[m/d]，dsim是简化裂缝段的开口大小，h1和h2是两端节点的水力头，Lsim是简化裂缝段的长度。对于图6a和6b所示的情况，n=1；对于图6c，n等于合并裂缝的数量。K通过裂缝水力传导率与开口大小之间的经验关系来估算：(2)K=AdsimB对于理想的平面裂缝，B=2，A=9.46。在本研究中，根据研究区域内裂缝的钻孔水力测试结果（刘等人，2024），采用经验值B=2.39和A=2.48。结合（1）和（2）得到：(3)dsim=Lsim?∑i=1nQori,iAh1?h2B+1为了评估基于流动约束的简化方法的优势，本研究引入了一种不需要流动约束的参考方法。在这种方法中，根据水力阻力来校正裂缝开口，以在没有流量数据（流量和方向）的情况下更好地保留原始裂缝网络的渗透能力。水力阻力R与裂缝的渗透能力成反比（De Hoop等人，2022）：(4)R=LKd=LAdB+1对于图6a和6b所示的情况，渗透能力得到保持（Rsim=∑Rori），简化段的开口大小计算如下：(5)dsim=Lsim∑i=1nLori,idori,iB+1B+1其中dori和Lori分别是原始裂缝段的开口大小和长度。对于图6c所示的情况，合并裂缝的等效阻力表示为：(6)1Rsim=∑i=1n1Rori,i将方程（4）代入方程（6）得到等效开口大小为：(7)dsim=Lsim∑i=1ndori,iB+1Lori,iB+1这种基于阻力的开口校正提供了一种不依赖于流场信息的直接方法来保持网络的整体渗透能力。它可以在基于几何形状的简化或预处理步骤中高效实施，确保简化后的网络具有水力上的一致性。2.3. 离散裂缝-基质模型数值建模的目标不是在所有局部尺度上都实现精确的数值再现，而是在可接受的误差范围内捕捉地下水流动和溶质传输的关键宏观水力响应。为了在更真实的条件下评估简化网络的稳健性，在验证阶段引入了一个离散裂缝-基质模型。在这个框架中，将简化的裂缝骨架嵌入到完整的裂缝-基质系统中，使我们能够量化裂缝-基质交换对网络简化引起的预测误差的影响。重要的是，这一步并不修改简化策略本身，而是用于评估其在流动不是严格由裂缝主导的系统中的适用性和局限性。这种方法阐明了算法的潜在限制，并确保简化网络在更复杂的水力条件下仍然可靠。比较了原始模型和简化模型之间的几个关键指标，包括总出流量、局部流场分布、溶质浓度空间演变和溶质突破曲线（BTCs）。这些指标的定量比较允许直接估计由简化引起的误差，并评估所提出方法在更真实的水文地质条件下的可靠性。在这个模型中，裂缝嵌入在多孔基质中。连续性条件确保压力和质量通量在裂缝-基质界面处一致。裂缝与相邻基质之间的正常方向上的通量交换被作为控制方程中的源项（Lin等人，2016）。研究域使用非结构化三角有限元网格进行离散化（图7b）。下载：下载高分辨率图像（254KB）下载：下载全尺寸图像图7. (a) 边界条件，以及 (b) 建模域的网格离散化。基质和裂缝中的流体流动可以分别表示为：(8)Matrix?um=0um=?kmμ?pm(9)Fractures?uf=Qfdfuf=?kfμ?pf其中um和uf分别是基质和裂缝中的速度向量（m/s）；km和kf分别是基质和裂缝的渗透率（m2）；pm和pf分别是基质和裂缝中的压力（Pa）；Qf表示裂缝-基质交换通量（m/s）。这里，流体是水，密度ρ=1000kg/m3，动态粘度μ=1×10?3Pa·s；2g=9.8m/s2是重力加速度向量。此处求解的流场作为直接溶质传输模拟的输入。对流-扩散方程描述了溶质通过基质和裂缝的传输过程：
(10) 基质：?εmCm?t?D?2Cm+um?Cm=0
(11) 裂缝：?εfCf?t?D?2Cf+uf?Cf=Jfdf
其中Cm和Cf分别是基质和裂缝中的浓度；εm和εf分别是基质和裂缝的孔隙度；D是有效扩散-弥散系数（m2/s）；Jf表示从裂缝到周围基质的溶质质量交换通量（kg·m?2·s?1）。

对于裂隙地下系统的数值模拟，在左侧和右侧施加了零通量边界条件，并在整个系统上施加了规定的压降，从而产生了从上到下的水力梯度（图7a）。此处求解的流场是溶质传输模拟的基础。初始和边界条件表示如下：
(12) CΩ=0 t=0
(13) CSin=C0=1 t≥0
(14) ?CSout?n=0 t≥0
其中n表示出口边界的法线方向。入口表面是一个狄利克雷边界，假设阶跃注入，C0=1.0（kg/m3），而出口表面是一个开放边界，没有扩散通量。

出口表面的通量加权 breakthrough 曲线（BTCs）是通过流出溶质质量与流体质量的比率计算得出的：
(15) Cf=∫?1uCdz∫?1udz

为了便于比较结果，浓度和时间进行了归一化处理：
(16) C′=Cf/C0
(17) t′=u?t/Lm
其中C′是归一化浓度，t′是无量纲时间或孔隙体积（PV），?u?是裂缝的平均流速（m/s），Lm是岩石基质的长度（m）。

均方根误差（RMSE）量化了简化模型和原始模型 BTCs 之间的差异：
(18) RMSE=∑?N(Csim′?Cori′)2
其中N是 BTCs 中的数据点数量，Csim′和 Cori′分别是简化模型和原始模型的归一化通量加权浓度。

流动和溶质传输的控制方程使用有限元软件 COMSOL Multiphysics 进行求解，该软件已广泛应用于裂隙多孔介质的数值模拟（Favier等人，2024；Li和Li，2019；Sun等人，2021；Xue等人，2023）。具体来说，使用了Darcy’s Law（dl）接口来计算流场，以及Transport of Diluted Species in Porous Media（tds）接口来模拟溶质传输。线性方程组使用直接求解器‘PARDISO’进行求解。所有模拟都在台式工作站上进行（Intel Core i7–11700 CPU，32 GB RAM，Windows 10）。在裂缝附近应用了精细的网格，最大元素大小为10 mm，以解决裂缝内及其周围的流动和浓度场的陡峭梯度，从而确保模拟的准确性。

2.4. 建模方案和参数设置
模型的计算效率取决于网格质量。在相同的网格设置下，更小的几何特征需要更细的网格。采用裂缝边缘长度阈值（lv）作为网络简化的控制参数。参数lv直接衡量了应用于网络的风裂尺度过滤。lv的值为20、40、60、80、100和120 mm，范围是根据原始网络中裂缝段长度的统计分布确定的。20 mm的下限略大于每个网络中最小的裂缝边缘长度。在这个阈值下，长度小于20 mm的裂缝段分别在网络A、B和C中占总段数的0.05、0.15和0.25。120 mm的上限是三个网络中最大的平均裂缝边缘长度，对应的比例分别为0.54、0.84和0.90（表1和补充图2）。相同lv在不同网络中产生的覆盖比率变化反映了裂缝长度分布和结构复杂性的内在差异。结构更复杂的网络包含更大比例的短裂缝段。

建立了一种无流动约束的简化方法作为基准，以评估所提出的流动约束的贡献。在遵循相同的几何简化程序（包括裂缝平滑、短裂缝删除和锐角消除）的情况下，基准方法有两个关键的方法学区别。首先，其简化决策省略了对流动方向影响的验证，仅依赖于最小化结构改变的几何或拓扑标准。其次，通过应用基于等效阻力的公式（(5)，(7)来实现水力补偿。这种实验设计使得能够隔离和定量评估流动约束在保持水力特性方面的贡献。

详细的参数设置总结在表2中。所有参数都是基于已发表的裂缝岩石实验选择的，使用现场测量值或代表性平均值（Geiger等人，2010；Wang和Cardenas，2014；Yuan-Hui和Gregory，1974）。

表2. 模型输入参数
参数 | 符号 | 值 | 单位
-------|-------|------|
| 水力梯度 | i | 1 | m/m |
| 动态粘度 | μ | 1×10?3 | Pa·s |
| 密度 | ρw | 1000 | kg/m3 |
| 注入浓度 | C0 | 1.0 | kg/m3 |
| 裂缝渗透率 | kf | 变量 | m2 |
| 裂缝孔隙度 | θf | 1.0-F | |
| 裂缝扩散率 | Df | 2.03×10?? | m2/s |
| 基质渗透率 | km | 10?1? | m2 |
| 基质孔隙度 | θm | 0.2 | |
| 基质扩散率 | Dm | 2.03×10?? | m2/s |

3. 结果和讨论
3.1. 简化网络几何
图8展示了在不同流动约束简化程度下裂缝网络拓扑结构的变化。相应的密度和拓扑参数列在表2中。与包含大量节点（红点）的原始网络相比，随着阈值lv的增加，节点密度和裂缝段数量显著减少。由于许多小而复杂的结构逐渐消失，剩下的主要骨架变得可识别且几何上更简单、更平滑。

3.2. 校正后的流场
在流动约束简化过程中，裂缝-基质系统内的流场使用(8)、(9)进行计算，如图9所示，对于不同的lv值。如图9a所示，简化改变了局部水力压力场（h），尤其是在移除或合并的裂缝交叉点和短段附近。这种几何上的减少自然减弱了细尺度的压力不均匀性，导致局部压力分布更加平滑。重要的是要注意，这种局部平滑是小尺度几何特征移除的直接结果，并不单独验证简化的有效性。在以裂缝为主导的系统中，小尺度压力不均匀性主要由几何不规则性（例如短段和复杂交叉点）控制，而不是由主导水力路径的变化控制。因此，细尺度压力波动的减弱并不一定意味着主导流动连通性或整体渗透性的丧失。关键标准是宏观上相关的水力行为——如主导流动结构和整体流动响应——是否得到保持。

3.3. 流动约束简化对溶质传输的影响
上述结果表明，流动约束简化可以在很大程度上保持原始裂缝网络的水力特性。然而，这种流场的一致性并不一定确保溶质传输行为的准确再现，因为溶质迁移不仅受主要裂缝渗透率的影响，还受局部速度、水力连通性和溶质扩散性的影响。为了研究网络简化对溶质传输的影响，基于稳态流场和流动约束简化算法模拟了不同简化程度（lv）下的溶质浓度和相应的 breakthrough 曲线。所有模拟在溶质浓度达到稳态后终止，并对简化和原始模型应用了相同的时间步长。

图12显示了原始和简化裂缝网络中的溶质浓度（t’ = 0.15PV）。与图9所示的主要流动通道相比，裂缝网络中的溶质传输在主要路径上表现出更多多样性。与允许高流量通过的主要流动通道不同，次要溶质路径表明在极低速度的裂缝中扩散占主导，特别是在垂直于主要流动方向的小裂缝中。简化程序清楚地消除了许多次要传输路径，保留了主要的骨架裂缝。随着裂纹密度的增加，主导通道的数量也会增加，例如从A到B和C的网络。然而，与图9中的lv相比，主导流动通道的变化很小，而主要的溶质传输通道随着lv的增加而减少，特别是对于更复杂的网络B和C。下载：下载高分辨率图像（537KB）下载：下载全尺寸图像图12. 使用流动限制简化方法（在0.15PV下）时，不同lv值的裂纹网络中的溶质浓度变化。图13a–c显示了出口表面（域底部）的标准化溶质通量（C′）的突破曲线（BTCs）。每个图中的子图1? C′详细说明了结果，特别是对浓度长尾的部分。图13d显示了简化网络和原始网络之间BTCs的RMSEs。总体而言，在所有情况下，BTCs都表现出一种典型的模式，即在突破后迅速上升，随后逐渐过渡到一个相对稳定的阶段。随着裂纹网络从A到B和C变得越来越密集，BTC的上升部分变得更加平缓，达到稳定阶段所需的标准化时间显著增加。下载：下载高分辨率图像（435KB）下载：下载全尺寸图像图13. 使用流动限制简化方法时，不同lv值的A、B和C网络的出口溶质通量的突破曲线（BTCs）以及RMSEs（d）。整个系统中的溶质传输是由裂纹中的对流和基质中的扩散共同控制的。裂纹中的流速通常比基质中的流速高1-5个数量级，这导致在传输的早期阶段，基质中的溶质浓度显著低于裂纹中的浓度。这种裂纹-基质界面处的浓度梯度不断驱动溶质从裂纹扩散到基质中。裂纹网络的简化主要影响两个物理过程：首先，溶质变得更加集中，并沿着少数主要通道传输，从而加快了向出口的整体迁移速度；其次，总的裂纹表面积减小，限制了溶质向基质中的扩散量。以裂纹网络C为例，与原始网络相比，简化网络的裂纹密度分别降低了75%、39%和27%（对于lv = 20、80和120）。随着简化程度的增加，流量逐渐集中在少数主要通道中，裂纹网络中每个节点的平均流速分别增加了1.3、2.17和3.01倍。最大流速分别增加了1.01、1.05和1.14倍，从而加快了沿主要通道向出口的溶质传输速度。在裂纹-基质质量交换方面，网络简化显著影响了溶质向基质的扩散。在标准化时间0.15PV时，原始网络中保留的总输入的比例分别为69.2%、69.0%、60.0%和53.0%，而lv = 20、80和120时，相应的累积溶质质量分别为10.0%、9.8%、7.6%和6.5%。这表明在简化网络中，大多数溶质在传输的早期阶段被保留在了基质中。随着传输的进行，基质中靠近裂纹的区域逐渐接近局部饱和，导致裂纹和基质之间的浓度梯度减小，进一步减少了溶质向基质的扩散。在1PV时，基质中累积溶质的比例分别降低到41.0%、40.3%、35.2%和30.0%，表明在简化网络中，更多的溶质通过裂纹传输到出口。在传输的后期（PV > 1），系统接近饱和，保留的溶质几乎可以忽略不计。由于简化网络的总体渗透性与原始模型的渗透性相差不到3%，因此此阶段的BTC与原始网络没有显著偏差。对三种裂纹网络结构的BTC曲线进行比较进一步验证了上述结论。具有简单裂纹结构和集中流的网络A，在早期阶段的出口浓度迅速上升。随着网络结构变得更加复杂（B和C），流动逐渐分布在更多路径上，溶质传输不再局限于少数主要通道，导致更多的溶质扩散到基质中。因此，随着简化水平的提高，模拟误差逐渐增加。尽管流动限制简化仍然高估了出口溶质通量的BTC，但如图14所示，这些误差比无约束简化方法的误差要小得多。具体来说，对于简单的网络A，当lv>80mm时，与原始模型的差异变得明显；而对于更复杂的网络B和C，当lv仅超过40mm时，与原始模型的偏离急剧增加。这些发现表明，如果没有水力场约束，过度去除或不适当的裂纹合并会改变主要通道上的流动路径和平衡，从而显著影响溶质传输的对流和扩散。下载：下载高分辨率图像（367KB）下载：下载全尺寸图像图14. 使用流动无约束简化方法时，不同lv值的A、B和C网络的出口溶质通量的突破曲线（BTCs）；RMSEs（d）。3.4. 简化裂纹网络的DFM的计算效率DFM的计算效率取决于划分的元素，例如三角形元素的数量（N???）和裂纹内的网格顶点数量（Nf），这些与本研究中使用的边长lv的网络简化指标高度相关。此外，还评估了两个其他指标：三角形元素面积比（R?）和平均元素质量（Save）。R?表示最小和最大元素面积之间的比率，反映了网格的均匀性。Save的范围是0到1，值越高表示元素形状更好，通过以下公式计算：(19)S=1?maxβ?βe180?βe,βe?ββ其中β是元素顶点的内部角度，βe是等边三角形中对应的角度（βe= 60°）。表4表明，边长lv的增加会减少三角形元素的数量（N???）和裂纹内的网格顶点数量（Nf），裂纹网络简化程序的时间消耗（Tsim）以及总模拟时间（T）。当lv从0增加到20mm时，网格元素的数量显著减少，因为原始裂纹网络（lv = 0）的网格分辨率主要受众多裂纹节点间距（0.1 ～ 2mm）的限制，而在本研究中则受到最大元素大小（10mm）的控制。同时，更高的简化程度（或更大的lv）通常会增加元素面积比（R?）和平均元素质量（Save）。因此，简化导致网格大小的分布更加均匀，接近理想的几何形状，最终提高了数值模型的收敛性和稳定性。随着简化水平的提高，简化程序的计算成本降低，因为更多的短裂纹段被直接移除，从而减少了拓扑检查和流动一致性迭代的次数。尽管对于结构复杂的网络，算法需要稍微更多的时间，但总体上的模拟效率提高显著。表4. 随着简化程度lv的变化，网格特性和模拟时间的变化。网络lv(mm)N???NfR?S???Tsim(s)T(s)A08018433190.0020.837-114720296436960.1300.8990.175940289276820.1240.9140.255560284676800.1240.9210.195480278436230.1320.9240.2254100277255950.1340.9250.2853120272835600.1210.9310.3152B0261958136820.0010.836-1901203017520840.0690.8784.90181402944217990.0530.9016.14162602876316660.0960.9026.98156802819313340.0140.9135.691401002726110950.1370.9245.301201202723110540.1070.9225.61118C0307982162160.0010.825-2292203061421540.0860.87713.14211402953518060.0630.90013.85179602833114440.0630.91213.16156802790211320.1060.91912.34139100269669110.1070.92511.69122120267827550.1350.93310.841144. 实地验证和方法局限性4.1. 实地验证：Chengqi露头案例本节的目的是验证所提出的裂纹网络简化方法在真实场地条件下的适用性和有效性。选择了一个尺寸为1m × 3.8m的完整暴露露头剖面（“Chengqi”剖面）作为实地案例研究。图15显示了露头的原始图像、提取的裂纹骨架以及相应的孔隙分布特性。裂纹网络的孔隙范围为0.4–2.4mm，裂纹密度为26m/m。岩石的岩性是中三叠纪关岭组的石灰岩（T?g）。数值模拟中使用了第2.4节中使用的相同物理参数，唯一的修改是将流动方向从左改为右。简化后的裂纹网络及其相应的头部分布结果在补充材料中提供。下载：下载高分辨率图像（710KB）下载：下载全尺寸图像图15. (a) Chengqi剖面露头和(b) 手动解释的裂纹网络。表5比较了简化前后网格大小、计算时间和溶质传输误差的变化。原始DFM模型需要1329,392个三角形网格元素进行离散化，溶质传输模拟时间为10,812秒。应用所提出的简化算法后，即使使用保守的参数设置（lv = 20mm），网格数量也可以减少到原始模型的2.8%，溶质传输模拟时间缩短到原来的5.4%，RMSE仅为0.008。这些结果表明，所提出的方法可以在保持高模拟精度的同时显著提高高分辨率裂纹模型的计算效率。表5. Chengqi模拟的数值性能。网络lv(mm)N???NfR?S???Tsim (s)T(s)RMSEChengqi01329392675560.0000.837-10812-203752757070.0200.8172125850.008403158041660.0280.8481653830.013603038132110.0810.8721223240.032802921827560.0330.881992850.0531002838723950.0300.888872440.0631202787719890.0840.897751990.0894.2. 简化阈值的选择简化阈值lv是一个关键的控制参数，它决定了几何减少和拓扑调整的程度。选择简化阈值lv受到三个主要考虑因素的指导：(1)拓扑稳定性和流动约束：裂纹去除和端点调整必须保持周围裂纹的主导流动方向，并保持物理上合理的全球拓扑——这是有效孔隙校正和质量守恒的先决条件。过大的lv可能会导致显著的几何变形，可能导致非物理性的重叠或交叉，从而损害网络的水力完整性。对所有测试网络的实证观察表明，lv的上限应使得修改后的裂纹累积长度不超过总裂纹长度的90%。这个限制确保了显著的计算收益，同时避免了对关键流动路径和宏观结构的不合理扰动。(2)算法性能：高分辨率裂纹网络模型的计算效率不仅受到拓扑复杂性的限制，还受到大量节点和小尺度几何特征的影响，这些特征在数值离散化中引入了过度的自由度。在低简化阈值（例如lv =20mm）下，算法主要减少节点密度和平滑裂纹几何形状，有效地消除了小的和复杂的局部特征，同时在数值效率上获得了显著的提升，而对整体网络拓扑的干扰最小。随着lv的进一步增加，并逐渐将更多裂纹段纳入简化过程，计算效率继续提高，但由于网格密度主要由全局离散化设置控制，而不是细尺度几何特征，因此收益趋于减少。同时，简化算法本身的运行时间也减少了，因为已经消除了许多短裂纹后，所需的拓扑检查次数减少。(3)对溶质传输的影响：在低简化阈值下，拓扑修改最小，裂纹网络的总体连通性基本保持不变。因此，溶质传输的准确性仍然很高，表明对微小几何特征的调整对宏观传输行为的影响有限。随着lv的增加，次级流动路径逐渐被移除，裂纹-基质界面面积减小，导致在主导主干裂纹内的对流更加集中，与基质的扩散交换减少（图13）。这种结构转变可能导致突破稍早发生，以及后期尾部的减弱。对于这里研究的网络，一个适中的阈值（大约lv = 40mm）达到了一个良好的平衡，在保持合理的RMSE范围内（RMSE < 0.02）的同时，显著降低了计算成本。总之，lv的选择应与预期的建模目的相匹配。当需要准确保留主导水力结构和传输行为时，相对较低的阈值（例如20–40mm）可以有效去除几何冗余，同时保持流动路径的稳定性。对于强调计算效率或大规模场景分析的应用，可以采用更高的lv来进一步减少结构复杂性，只要在可接受的误差容忍范围内保留了主导的水力响应。4.3.方法的局限性：来自四个断裂网络的模拟结果表明，所提出的简化方法能够在测试范围内保持稳定的水力响应，无论网络规模或结构复杂性如何。该算法通过局部几何平滑和拓扑更新来运作，而不是依赖于整个网络的大小或节点密度。因此，即使在不同的空间尺度下，也能一致地保留主要的控流路径。然而，需要注意的是，尽管相对效率有所提高，但绝对计算成本仍然会随着模型规模和断裂密度的增加而增加。该方法基于一个单一的稳态流场构建，其有效性取决于主要流动模式保持相对稳定的假设。当边界条件发生显著变化时（例如，季节性补给变化或泵送引起的梯度反转），优先流动路径可能会重新组织。在这种情况下，过高的简化阈值可能会使网络过度依赖于初始的水力状态，从而降低其在其他流动场景下的鲁棒性。因此，对于流动方向或梯度幅度可能发生变化的应用，建议使用保守的阈值（例如lv ≤ 20 mm）。这种方法将结构修改限制在较小的几何调整范围内，从而在保持主要传导框架的同时，仍提供适度的计算改进。

本研究中的所有案例研究都基于二维（2D）断裂网络。实际上，详细的断裂几何信息主要来自2D来源，如露头调查、地表测绘和井洞成像，而能够完全约束三维（3D）断裂连通性的野外数据集仍然很少。因此，开发适用于高分辨率2D断裂网络的稳健简化方法具有直接的工程意义。所提出方法的基础图论概念原则上可以应用于三维离散断裂网络，在这种网络中，断裂可以被表示为多边形元素，其顶点和交点可以被抽象为图节点和边。然而，将算法扩展到3D需要在几何处理、交点检测和数据结构管理方面进行大量修改，这仍然是未来研究的课题。

5. 结论：从喀斯特露头剖面中提取的高分辨率断裂网络为地下水流动和溶质传输的建模提供了详细的几何和拓扑表示。然而，这些网络的结构复杂性往往导致庞大的数值模型和高计算成本，限制了基于场景的模拟和不确定性分析的可行性。为了解决这一挑战，本研究提出了一种基于流动约束的断裂网络简化方法，并系统地评估了该方法在不同规模和结构复杂性网络上的性能。结果表明，所提出的方法在保持主要流动路径和整体水力响应的同时，大幅降低了网格密度。在不同简化水平下，等效渗透率保持在狭窄的偏差范围内，突破曲线差异也保持在可接受的限度内。与无约束的几何简化方法相比，基于流动约束的策略更有效地保持了关键的传导结构和传输行为。进一步建议根据建模目标采用分层阈值策略：较高的阈值适用于计算密集型应用，可以快速提取主要的传导骨架；而较低的阈值则适用于强调详细溶质传输或断裂-基质相互作用的研究。在这个框架中，简化阈值的选择应基于水力性能标准，而不仅仅是几何考虑。

提出的简化方法基于对断裂网络的几何和拓扑调整，其标准不依赖于特定的空间维度，因此可以扩展到三维断裂网络。尽管该方法在当前的单相饱和流动条件下已经显示出显著的计算加速效果，但其是否适用于更复杂的物理过程（如多相流动和反应性溶质传输）仍需要进一步验证。总体而言，所提出的方法为量化并表征喀斯特裂隙储层中的流动和溶质传输不确定性提供了一个高效、稳健且可扩展的框架。

作者贡献声明：
陈曦（Xi Chen）：撰写原稿、方法论、资金获取、形式分析、概念化。
高凤军（Fengjun Gao）：撰写原稿、可视化、软件、方法论、调查、形式分析、数据管理。
张志才（Zhicai Zhang）：撰写原稿、方法论、资金获取。
刘伟汉（Weihan Liu）：资源支持、概念化。
梁英春（Yingchun Liang）：撰写原稿、方法论。