岩溶管道网络数值模拟方法研究进展

王益伟，渠光华，刘 埔，褚学伟，段先前

（1.贵州大学 资源与环境工程学院，贵阳550025；2.贵州大学 矿业工程学院，贵阳550025）

岩溶管道网络控制着岩溶含水层中地下水流的水力特征及运移特征，对岩溶含水层中地下水流及溶质运移过程有着重要的影响［1-2］。要清晰刻画地面以下岩溶管道网络，需要获得大量精确的管道网络的空间信息［3］，然而，由于管道网络的复杂性及非均质性，导致实际工作中是无法直接获得充足的空间信息的。因此，采用数值模拟的方式获得缺失的空间信息，是重现真实管道网络的有效手段［4］。

目前，岩溶管道网络特征的数值模拟研究，主要集中在：①岩溶管道网络系统演化过程模拟；②岩溶管道网络结构模拟（Structure-imitating）。

1 岩溶管道网络系统演化过程模拟

岩溶管道网络系统演化模拟，其目的是模拟控制碳酸盐岩溶解的完整物理化学过程，包括模拟空隙、管道分布特征与水流及可溶岩溶解过程的关系，最终实现对整个岩溶管道系统演化的物理和化学过程的模拟［5-6］。这种方法对于清晰的理解岩溶演化过程是非常有用的，并能够提供较为精确的关于空隙、管道的演化信息。但该种方法较难应用于实际中，很难反映出真实管道网络系统的空间特征，原因在于重建形成当前管道系统的古气候及地质条件是相当困难的。同时，构建这样的模型，需要大量的数据，并且该类模型的解算过程通常是非线性的，对CPU的要求相当高。

2 岩溶管道网络结构模拟

岩溶管道网络结构模拟研究是一个极具挑战的任务。原因在于岩溶管道网络的复杂性，致使客观上不可能存在有足够详细的信息去详尽完整地反映出实际局部岩溶管道的特征［7-8］。针对这一问题，一种思路认为解决真实局部岩溶含水层的水文地质特征的预测，不如去解决理解整个岩溶含水系统的功能问题［9］。另一种思路，则采用统计学的观点来重现管道网络［10］，从而延伸弥补缺乏空间信息的区域。

2.1 全局式模型

基于第一种思路，全局式的模型被提出。全局式的模型，以黑箱理论为基础，将岩溶管道网络系统作为一个整体，利用系统的全局式反映，如利用泉的水动力及水化学特征，推测反演岩溶管道网络特征。如，根据泉点流量系列、水化学系列及同位素数据，使用时间序列分析、主成分析等方法，建立起多种解译模型［11-12］。Mayaud等［13］通过对泉点的时间序列分析，建立了有限元模拟模型，阐述了管道网络对研究区域的水文地质条件的影响。长勇和吴吉春等［2］以泉流量系列为基础，利用MODFLOW-CFP模块构建了岩溶含水层的管道网络模型，并分析了管道网络对泉的流量过程的影响。尽管这些方法被广泛的采用，但全局式的模型只揭示了系统全局的信息，大多数情况下并不能清晰地描述岩溶管道网络的空间特征。这种模型能够解决区域水量及流向问题，但对小尺度下局部水流的流量及流向问题是无法刻画的，这样的模型实际上是无法满足需要明确知道水流运移路径的研究，如岩溶含水层污染物运移及修复的研究。

2.2 分布式模型

基于第二种思想，一些分布式模型被提出。分布式模型能够提供更真实的岩溶含水层中地下水的流量及水流方向的信息，但其难点及关键在于如何处理岩溶管道网络的非均质性，既在空间中的变化。方法包括：地质统计的方法、分形的方法［14］、格子气元胞自动机（lattice-gas cellular automaton，LGCA）方法［16］、伪并行遗传算法［2］及反演模拟［3］。其中，基于统计学原理的模拟方法被广泛使用。其余几种方法通常会生成一个二维模型，但实际的管道网络应该是一个三维的空间实体，同时，在这些模型中一些影响管道网络形成的常规原则很难约束到所建的模型中。

基于地质统计学原理的模拟方法根据地质统计学原理，建立管道网络的空间几何特征（方向、相互连接方式、连通性等）与可野外直接测量得出的参数（裂隙、基底面、水力梯度、水流出口位置）之间的相关关系，从而实现对管道网络空间分布特征的模拟。其中，最常用的方法是序贯指示模拟，用指示随机变量来描述管道的空间变化［16-17］。Jaquet及Jeannin［15］使用管道的密度作为相关空间变量，采用转向带法法模拟了在某一区域上的管道网络，并根据模拟的密度值随机的确定了管道的位置。但以上这些传统模型都基于变差函数的两点地质统计学方法。由于两点统计测度只描述在已知一点的值的情况下，在特定位置上特定值的可能性，不受两点之间结构及通过复杂的路径是否连通的影响，因此传统的地质统计模型不足以表达具有复杂空间结构的含水介质，是没有能力描述含水介质的连通性的［4，18］然而，对复杂的岩溶管道网络来说，连通特性无疑是一个非常重要的特性，是衡量岩溶含水层非均质性的一个重要参数。这种特性使得含水介质即便具有相同的网络结构，当其连通性不同时，水流及溶质在管道网络中的运移也会发生变化。为了解决模型中含水介质连通性问题，基于多点地质统计的高维度的非均质模型被提出。这些模型的基本思想是利用空间多个点的组合模式来描述含水介质结构信息，采用训练图像代替变差函数来描述地质体结构信息，克服两点地质统计不能再现目标几何形态的不足，这类型的模型更适合于进行具有复杂结构含水介质信息的模拟。采用多点地质统计方法的关键在于构建合理的训练图像，尤其是在模拟三维的地质体时［19-21］。训练图像构建的方法包括基于对象、 基于过程的方法及混合方法［20-22］。尽管，多点地质统计适用于模拟复杂的地质结构，同时岩溶含水层具有复杂的管道网络，但目前多点地质统计学方法多用在模拟裂隙介质及多孔介质［23-25］，在油田工程方面应用较多，却很少用于对岩溶管道网络的模拟。

2.3 非均质模型

伴随着复杂的非均质模型的出现，连通性的测度成为一个量化非均质模拟及标度非均质化模型的重要工具。连通性测度的研究主要基于两种理论：逾渗理论［26］及图论［27］。

以逾渗理论为基础所建立的连通性，需要事先知道渗流的关键阈值，取一个简单的渗透系数作为连通性的判断。但Edery等［28］的研究结果却表明低渗透地区也可能出现较大的流量，而一些低流速的区域并不连通。基于图论的连通性测度研究，对连通性的定义更为清晰。Renard等［4］和Tyukhova等［29］基于图论的原理提出了一个简单的定量网络连通性的方法，通过确定最小路阻的路径，最终确定裂隙含水介质的连通性。同时，在Renard等［4］的研究中将含水介质的连通性分为静态连通性及动态连通性两种。静态连通性取决于含水介质中存在高渗透性的通路，依赖于介质本身的渗透性能。而动态连通性则取决于水流及物质运移的物理化学过程，如水流及溶质的沟槽流现象。由于静态连通性不易被直接测出，但动态连通性较易测出，通常可通过野外的示踪实验获得，因此，揭示静态连通性和动态连通性之间的关系，有助于利用野外现场实验提高数值模型对真实岩溶管道网络的反演及再现，是管道网络模拟研究中的一个重要的内容。静态连通性和动态连通性之间存在复杂的关系［20-32］。如，Ronayne等［33］提出，当在岩溶管道周围出现高渗透性的围岩时，静态连通性不能保证存在高的动态连通性。Tyukhova and Willmann等［34］根据管道的路阻及管道宽度定义了静态连通性的测度，并将其与动态连通性进行对比，最终将其测度值应用于对水流及溶质运移的模拟。尽管，目前对岩溶管道网络连通性的研究不多，但该特性的研究将有助于提高模型对非均质含水层中水流及溶质运移模拟的预测能力；同时，量化岩溶管道网络连通性的测度，分析动态连通性与静态连通性的关系，是将衡量岩溶含水层非均质特性引入岩溶含水层水流及物质运移模拟研究关键方法。

3 结语

（1）岩溶管道网络模拟研究历经几十年的努力，取得了许多成果。早期研究多集中建立全局式、分布式的模型，忽略含水层内部的复杂性，随着计算机技术及物探技术的不断提高，要求更高精度的非均质模型将成为解决岩溶含水层由于其复杂性所导致的工程问题的有用工具。

（2）在非均质模型中，目前多基于地质统计学方法构建管道网络模型，可更真实地描述岩溶含水层中地下水的流量及水流方向的信息，但传统的地质统计模型多基于两点地质统计学方法，其所构建的岩溶管道模型不足以描述岩溶含水层中管道空间结构的复杂性，缺乏描述具有强烈非均质性含水介质连通性的能力。采用基于多点地质统计方法，结合空间中多点特征，综合反应含水介质空间中的结构信息，构建能够反映管道网络连通特征的管道网络模型，将是管道网络数值模拟的新趋势，其将对解决预测岩溶含水层中优势水流路径及溶质晕体分布特征等关键问题有着重要意义。