岩溶地下水系统演化的数值模拟

王 喆，卢 丽，夏日元，石胜良

(1.中国地质科学院a.岩溶地质研究所;b.国土资源部岩溶动力学重点实验室，广西桂林 541004;2.广东省地质矿产公司，广州 510083)

1 研究背景

我国岩溶地区发育面积较大，约占全国面积的1/3左右，在这些地区发育着各种泉，这些泉在我国经济发展的过程中，发挥着重要的作用。因此需要准确地模拟出泉流量，这样才能为泉域的水资源地保护和开发利用提供必要的依据。目前，建立泉流量模型一般有4种方法:线性法、非线性法、数值法和耦合法。线性法一般包括多元回归法［1］、偏最小二乘法等;非线性法包括灰色预测模型法［2］、神经网络法［3］、时间序列分析法［4］等;数值法主要是运用MODFLOW软件进行预测分析［5］;耦合法则是目前较为流行的方法，有的是线性法和非线性法的耦合［6］，有的则是非线性法和数值法的耦合［7］。从上述论述可以看出，耦合法可以将2种或者多种方法进行融合，从而克服因算法本身缺陷所产生的问题，可以提高模型的精度，因此耦合法是泉流量模拟研究的发展趋势。目前的耦合模型较少关注泉流量与岩溶含水系统演化的关系，只是追求泉流量的拟合效果，而且从未有过学者在分析岩溶含水系统的演化与泉流量关系中考虑过岩溶溶蚀过程，但溶蚀过程却是一个重要阶段，不能忽略。

因此本次研究利用基于VC++6.0的Fracture-ToKarst程序对泉域系统的渗流模型和溶蚀模型进行耦合，分均质和非均质2种情况构建模型，并对泉流量进行分析。此程序适用于非连续裂隙网络介质模型，具体耦合过程是首先利用已知条件对渗流场进行计算，得到每一段裂隙水头和流速，计算出钙离子的浓度，依据公式算出每一段裂隙扩展速度及一个步长后裂隙的新宽度，在新裂隙宽度条件下重新计算渗流场，并依据上述步骤进行下一个步长的裂隙宽度计算，从而使渗流模型和溶蚀模型达到耦合。

2 泉域系统演化数学模型构建

在构建耦合模型前必须选择渗流模型。目前，岩溶地区渗流模型有3种:连续介质模型、不连续介质模型和耦合介质模型。本文选用非连续裂隙网络介质模型。该模型认为裂隙与岩体间不存在水量交换，水流只在等效裂隙网络中运动［8］。

模型的构建过程大致如下:首先将研究区的裂隙分为确定性裂隙和随机裂隙2种，确定性裂隙直接输入模型，随机性裂隙通过Monte-Carlo产生，这2种裂隙构成原始裂隙网络。将没有导水能力的裂隙进行剔除，构建连通型裂隙网络，接下来对裂隙网络进行离散化处理。裂隙网络建好后，利用渗流模型计算出裂隙节点水头，构建初始的渗流场，然后利用初始渗流场和裂隙网络来计算岩溶系统中裂隙的溶蚀过程，推演泉域系统演化过程，进而预测泉流量。

2.1 模型假设

(1)介质场为等效裂隙网络模型，由确定性裂隙和随机裂隙构成，确定性裂隙直接输入，随机裂隙由Monte-Carlo法生成，裂隙两两正交，裂隙段假设为光滑平板裂隙。

(2)渗流场中的水流为稳定流态。

(3)溶蚀场的碳酸盐溶蚀速率，利用White的半经验公式进行计算。

2.2 裂隙网络模拟

裂隙网络分确定性裂隙和随机裂隙。对于确定性裂隙，根据长度、宽度、产状等参数，构建确定性裂隙网络，对于数量众多的随机裂隙，采用 Monte-Carlo模拟技术［9］生成随机裂隙网络，将确定性裂隙网络和随机裂隙网络耦合在一起，统一编号，然后对统一后的裂隙网络进行离散化处理。

2.3 裂隙网络中渗流场模拟

根据水流为稳定流态和裂隙段为光滑平板裂隙这2个假设条件，可得出其水流方程满足立方定律，即

式中:q是单宽流量;ρ是质量密度;μ是流动黏滞系数;b为裂隙张开度;H是裂隙两端节点的水头。

按水量平衡原理，节点i处的稳定流水流方程为

式中:qj为裂隙j单元流入(为正)或流出(为负)节点i的流量;Qi为节点i处的源汇项。

将式(1)和式(2)联立起来，采用迭代方法进行求解，进而得出渗流场。

2.4 岩溶裂隙网络溶蚀模拟

由模型假设可知，本次研究裂隙形态是基于平行板模型的单裂隙，因此本文采用White［10］所提出的单裂隙溶蚀扩展数学方程，即

式中:F是溶蚀速率(mol·cm-2·s-1);Kc是溶蚀速度常数(cm4·mol-1·s-1);Ceq是Ca2+的饱和浓度;C是水中某时刻Ca2+浓度。

在一段裂隙中，溶液带走的钙离子的数量可以表示为

式中:Q是每段裂隙中的流量;Δt是溶蚀作用的时间步长;C0是溶液流入一段裂隙时的初始钙离子浓度。

由质量守恒定律和式(4)可以推算出，单裂隙张开度扩张速度为

式中:ΔB是每一个步长内裂隙张开度变化量;M是在一段裂隙中溶液带走的钙离子的数量;ρ1是可溶岩密度;L是裂隙段长度。

如图1所示，对于一个裂隙均衡单元，假如裂隙交叉点i0的水头及Ca2+浓度为H和C，其相邻交叉点 in(n=1，2，3，4)的水头及 Ca2+浓度为 Hn，Cn，则裂隙交叉点i0与其邻点in有

其中:

式中:Ci为邻点in的钙离子浓度;ΔCi为水流从节点in到节点i0的钙离子浓度增加量;qi为节点in到节点i0的流量;Li为节点in到节点i0的裂隙迹长;Bi为节点in和节点i0的宽度。

上述公式的推导过程参见袁道先所写的《中国岩溶动力系统》［11］。对式(5)采用迭代法进行求解。具体步骤如下:首先利用式(6)至式(9)迭代算出节点Ca2+浓度，然后在根据式(4)和式(5)计算出裂隙张开度扩张速度。

图1 裂隙网络离散单元Fig.1 A unit in discrete fracture networks

3 泉域系统演化模拟与分析

在已有数学模型的基础上，结合补径排、边界条件和裂隙网络等情况，对泉域的演化进行模拟和分析。

补径排条件:假设模型只接受降雨入渗补给和河流侧向补给。这其中，降雨入渗补给的范围是图2(a)绿线部分，年降雨量为400 mm，降雨入渗系数为0.15，模拟期内降雨入渗补给量不变(本次模拟期为5 000 a);河流侧向补给的范围是图2(a)蓝线部分，水位高程为30 m，排泄河流高程为10.5 m，河流侧向补给量在模拟期内不变，且小于降雨入渗补给量。

图2 岩溶泉域系统概念模型Fig.2 Conceptual model of the karst spring system

边界条件:补给河流设定为定水头边界，补给河流与排泄河流之间的地面线处的边界设定为降雨入渗补给边界，泉水点设定为自由排泄边界，泉水点之间的节点也按自由排泄边界进行处理，而系统左边、底边及右下边界为隔水边界(如图2(a)中红线所示)。

裂隙网络:由确定性裂隙和随机裂隙组成的(图2(b))，图2中所显示裂隙网络的只是形象化显示，并非实际情况。真实的裂隙网络分均质系统和非均质系统2种情况，均质系统裂隙网络无确定性裂隙，由随机裂隙组成，且根据模型假设条件，随机裂隙两两垂直，且水平裂隙和垂直裂隙迹长分别为400 cm和200 cm，间距分别为200 cm和100 cm，表示为Lh和Lv，横纵向的张开度均为0.02 cm，裂隙条数分别为841条和735条;非均质系统裂隙网络由一条确定性裂隙和随机裂隙组成，其中确定性裂隙迹长为100 m，在侵蚀基准面处，张开度为0.08 cm，用ap来表示，而随机裂隙则与均质系统的相同。

3.1 均质性岩溶泉域系统演化

根据上节所介绍的裂隙网络情况，均质性岩溶泉域系统的裂隙网络是通过如下步骤进行构建:首先根据模型假设，均质系统的随机裂隙分为水平裂隙和垂直裂隙，然后根据上节所介绍的裂隙参数情况将水平裂隙和垂直裂隙两者的迹长、间距、倾角和张开度等参数输入FractureToKarst程序中，然后根据Monte-Carlo法构造各参数的随机数发生器，进而形成裂隙网络。

均质性岩溶泉域系统演化过程如图3所示，图内的红线为含水系统内等水位线。下面详述一下各个时间点上含水系统水位和张开度的变化情况。

图3 均质性岩溶泉域系统演化过程Fig.3 Evolution of a ho mogeneous karst spring system

100 a:岩溶泉域系统由于处于溶蚀的初期，其排泄区水头值较高，在右侧边界处形成了6个泉点排泄口，发育高程分别为 11，12，13，14，15，16 m，排泄河流高程为10.5 m，泉流量为0.37 mL/s。而在张开度方面，只在河流补给处，表层裂隙部分被溶蚀，张开度变大，深层裂隙没有大的变化。

1 000 a:排泄区水位较100 a时下降明显，下降约4～8 m，泉点变成2个，发育高程分别为11 m和12 m，排泄河流高程为10.5 m，泉流量为1.96 ×104mL/s，中部水位下降1～2 m，其他地方水位没有变化。而在张开度方面，排泄区和中部裂隙进一步被溶蚀，张开度变大，仍在表层范围内，没有向深层发展。

2 000 a:较1 000 a的情况，排泄区水位下降1～2 m，系统中部水位大幅度下降，下降约10 m，泉点仍为2个，发育高程分别为11 m和12 m，排泄河流高程仍为10.5 m，泉流量为2.27 ×106mL/s。在张开度方面，溶蚀范围已经由表层发展到了深层，基准面以下的裂隙大部分被溶蚀。

5 000 a:与2 000 a相比，含水层中部水位稍有下降，泉点仍为2个，发育高程分别为11 m和12 m，排泄河流高程仍为10.5 m，泉流量为3.90 ×108mL/s。而在张开度方面，裂隙扩张的范围已经发展到河流补给区。

由上述过程可以看出，2 000 a以后含水系统侵蚀基准面附近的裂隙大部分被溶蚀，这个现象在5 000 a时更加明显，溶蚀范围已经扩展到河流补给区。由于2 000 a后泉点始终维持在2个，因此随着溶蚀过程的不断发展，逐渐形成了侵蚀基准面附近的2个大型岩溶管道。由于岩溶管道的存在，进一步的加大了整个系统的溶蚀进程。

为了考察裂隙张开度溶蚀扩展情况，分别做横剖面图和纵剖面图(如图4)。如图4(a)(高为11 m的横剖面)所示，100 a时裂隙张开度几乎没有变化，仍为0.02 cm;到了1 000 a时，裂隙张开度略有变化，最大值为1.62 cm;2 000 a时，张开度变化较大，水平方向30～100 m内，张开度增大明显，最大值为9.69 cm，而0～30 m内则没有变化;4 000 a和5 000 a张开度增幅情况与2 000 a时大致相同，大约为 8 cm/ka，最大值分别为25.84 cm和33.91 cm，张开度变化范围则进一步扩大，为14～100 m。从2 000～5 000 a，张开度的增幅近似恒定，对比图3中后3幅图可以发现，这段时间内，裂隙被溶蚀的范围很大，并且在侵蚀基准面处形成了2个岩溶管道，汇集了大部分的补给水量，由于补给量恒定，因此裂隙张开度的扩张速率也是相同的。

图4 裂隙张开度演化情况Fig.4 Evolution of the fracture aperture

如图4(b)(长60 m的纵剖面)所示，1 000 a前裂隙张开度的变化不是很大，而2 000 a后张开度的增幅趋于恒定，为8 cm/ka，但14～30 m范围内的裂隙张开度几乎没有变化，原因是2 000 a后，岩溶管道在侵蚀基准面处基本形成，使得这些裂隙都无法被进一步溶蚀变宽。

3.2 非均质性岩溶泉域系统演化

本文在均质性模型的基础上，在侵蚀基准面处加上一条张开度为0.08 cm的大裂隙，以此构建非均质性模型。

如图5所示，跟均质性系统相比，非均质系统在演化初期水位下降较快，100 a时已相当均质性系统2 000 a时的水平，但在演化中后期水位几乎没有变化，这是因为侵蚀基准面处有一条张开度为0.08 cm的大裂隙，加速了系统的演化，并在1 000 a时发育成岩溶管道，水流大部分从岩溶管道流走，形成了稳定的排泄通道，水位得以相对稳定。泉点则由初期的2个变为中后期侵蚀基准面处的1个，演化初期泉点的高程为11 m和12 m，中后期则为12 m，而排泄河流的高程始终为10.5 m。

图5 非均质性岩溶泉域系统演化过程Fig.5 Evolution of a heterogeneous karst spring system

由于大裂隙的存在，差异性溶蚀现象更加明显，相同时间2种系统裂隙被溶蚀范围和张开度拓宽速度都有很大的不同。比如，2 000 a时非均质系统的溶蚀范围已延伸到补给区，比均质系统宽近30 m，但在深度方面则相对较浅。从图6中我们就可以更直观的看出非均质系统中裂隙张开度的变化情况，1 000 a后张开度的增幅速度，与均质系统增幅一致，但相同时段内张开度的最值却比均质系统要大4～5 cm。

图6 x=60 m处裂隙张开度纵剖面图Fig.6 Evolution of the fracture aperture(width vs.aperture)at the 60m length

4 泉流量预测研究

图7(a)、图7(b)2个子图分别是均质性及非均质性岩溶泉域系统泉流量曲线拟合图。从子图7(a)中可以看出，50 a至800 a这段时间，泉流量一直＜9 mL/s，而在900 a时，泉流量有了较大的增长，为1 265.16 mL/s，1 000 a 时增加到了19 590.45 mL/s，1 000 a后泉流量更是以指数函数的形式在增加。

产生上述情况的原因是由于岩溶含水系统具有差异性溶蚀现象，而侵蚀基准面具有较强的势汇，其逐渐开始汇聚较多的补给水流(包括降雨补给和河流补给)，这些水流携带有较多的CO2，侵蚀能力强，而且此处水循环深度浅，源汇距离短，径流强烈，因此逐渐在侵蚀基准面处形成了大型岩溶管道，不断地袭夺整个含水系统中的水量，这导致了1 000 a后泉流量以指数递增。

而深层原因是在1 000 a时，含水系统发生穿透现象，裂隙溶蚀已不再局限在表层范围内，而是向深层发展，张开度迅速增大，有几条裂隙都是由表层直接穿透到侵蚀面，而且侵蚀基准面处已经形成了小型岩溶管道，降雨量及河流补给量沿着贯穿裂隙和岩溶管道迅速排出系统，此时的优势通道最小裂隙张开度为2 cm左右，前人的研究将此作为岩溶系统演化前期和后期的分界线［12］，这就是1 000 a前泉流量增长缓慢，而1 000 a后指数增长的原因。

图7 岩溶泉域系统泉流量曲线拟合Fig.7 Fitted curves of spring discharge in the karst spring system

根据补径排条件可知，降雨入渗量为泉域主要补给来源，随着含水系统不断的溶蚀演化，水位不断下降，这使得降雨补给深度变短，当发生穿透现象后，降雨补给量比河流侧向补给量更快沿着贯穿裂隙到达岩溶管道，然后迅速排出系统，而河流侧向补给水量和范围有限，且影响范围内贯穿裂隙相对较少，因此泉流量主要来自降雨入渗补给量。由于泉水流出后汇入排泄河流内，因此排泄河流流量与泉流量成正比，即随着系统不断的溶蚀演化，排泄河流流量也随之增大，但是并不是无限增大，由于补给量一定，其溶蚀能力也是有限的，当侵蚀基准面处的岩溶管道发育到一定程度必然会停止，排泄河流流量也会趋于稳定。

从对泉流量演化机理的分析可以看出，影响泉流量的因素有以下几点。①补给水量:排泄量的大小与补给量密切相关，在岩溶地区更是如此，而本文含水系统在1 000 a时，发生了穿透现象，补给水通过贯穿裂隙和岩溶管道迅速流动，以泉的形式进行排泄，因此其成为了影响泉流量大小的重要因素。②张开度:由立方定律可知，流量与张开度的三次方成正比，1 000 a后，受穿透现象的影响，岩溶管道开始袭夺含水系统中的水流，而本文的补给量是恒定的，因此张开度的大小直接影响泉流量的大小。

根据泉流量的演化机理、影响因素和泉流量随时间变化曲线，推算出均质性系统和非均质性系统泉流量公式:

均质

非均质

非均质性系统泉流量变化曲线形状与均质系统的大体一致，只是其泉流量增幅起始时间提前(由1 000 a提前到了600 a)。相比于均质性系统，影响非均质系统泉流量大小除了补给量和张开度两个原因外，还有一个重要原因就是非均质性，由于演化初期存在着大裂隙，差异性溶蚀现象明显，大型岩溶管道形成早，穿透现象也出现得早(600 a)，张开度拓宽速度加快，导致泉流量增幅时间提前，同时段非均质系统泉流量比均质性系统的要大。

上述式(10)和式(11)只适用于裂隙岩溶介质系统，由于本文是在假设条件下对泉流量进行拟合，因此上述公式只是在假设条件下才会出现。本文假设裂隙两两正交，相同条件下，裂隙段内水流可以更快流动，使泉流量达到最大值，如果裂隙斜交，会使水流流动变慢，泉流量也会减小，因此裂隙相交的角度对泉流量会有一定影响，图7所示的泉流量在2个系统内均是逐渐增大的，但是并不会一直增大下去，本次模拟期为5 000 a，由于补给量是一定的，当岩溶溶蚀演化到一定阶段时，裂隙的张开度就会不在变宽，泉流量也会趋于稳定不在增长。

5 结论

本文在耦合渗流模型和溶蚀模型的基础上，对岩溶泉域系统进行了概化，分别模拟了均质性系统和非均质性系统的演化过程。均质性系统随着时间的推移，其排泄区和中部水位逐渐下降，裂隙溶蚀范围由表层发展到了深层，在2 000 a后在侵蚀基准面附近形成了2条大型岩溶管道;非均质系统与其发展模式类似，只是差异性溶蚀更加明显，在1 000 a左右的时间在侵蚀基准面处形成了大型岩溶管道。张开度的增幅速度2个系统大致相同，为8 cm/ka，但是起始时间不同，非均质系统比均质系统早了1 000 a，而且相同时段的张开度的最大值也要大4～5 cm。

随着时间的推移，2个含水系统不断被溶蚀，泉流量也在不断地增大。由于穿透现象和岩溶管道的出现，使得均质性系统1 000 a后泉流量有了大幅度的增加，而非均质性系统泉流量拐点出现在600 a。根据泉流量的数据对2个系统的泉流量曲线进行了相应的预测，均质系统为指数函数，非均质系统为多项式函数。

［1］姜宝良，付北锋，赵延涛.泉水动态分析预测和资源评价-以辉县百泉为例［J］.水文地质工程地质，2002，(3):45-48.(JIANG Bao-liang，FU Bei-feng，ZHAO Yan-tao.Analysis Prediction and Resources Evaluation of Spring Dynamic:A Case Study of Bai Spring in Hui County［J］.Hydrogeology and Engineering Geology，2002，(3):45-48.(in Chinese))

［2］郝永红，黄登宇，张文忠，等.山西神头泉流量的灰色预测模型研究［J］.水利学报，2004，(4):111-114.(HAO Yong-hong，HUANG Deng-yu，ZHANG Wenzhong，et al.Gray Prediction Model for Forecasting Spring Discharge［J］.Journal of Hydraulic Engineering，2004，(4):111-114.(in Chinese))

［3］陈宏峰，朱明秋，夏日元，等.湖南洛塔干河猪场表层岩溶泉BP人工神经网络分析［J］.中国岩溶，2005，24(4):300-304.(CHEN Hong-feng，ZHU Ming-qiu，XIA Ri-yuan，et al.Analysis on Epikarst Spring with BP ANN in Luota，Hunan Province［J］.Carsologica Sinica，2005，24(4):300-304.(in Chinese))

［4］李国敏，杨 福.神头泉流量衰减的时间序列分析［J］.中国岩溶，1993，12(2):123-130.(LI Guo-min，YANG Fu.Time Series Analysis for the Declining Discharge of Shentou Spring［J］.Carsologica Sinica，1993，12(2):123-130.(in Chinese))

［5］王庆兵，段秀铭，高赞东，等.济南岩溶泉域地下水流模拟［J］.水文地质工程地质，2009，(5):53-60.(WANG Qing-bing，DUAN Xiu-ming，GAO Zan-dong，et al.Groundwater Flow Modelling in the Ji'nan Karst Spring Area［J］.Hydrogeology and Engineering Geology，2009，(5):53-60.(in Chinese))

［6］陈南祥，黄 强，曹连海.基于偏最小二乘回归与神经网络耦合的岩溶泉预报模型［J］.水利学报，2004，(5):68-72.(CHEN Nan-xiang，HUANG Qiang，CAO Lian-hai.Model for Prediction of Karst Spring Flow Based on the Coupling of Neural Network Model with Partial Least Square Method［J］.Journal of Hydraulic Engineering，2004，(5):68-72.(in Chinese))

［7］鞠 琴，郝振纯，余钟波，等.泉流量模拟研究——以小南海泉域为例［J］.四川大学学报(工程科学版)，2011，43(1):77-82.(JU Qin，HAO Zhen-chun，YU Zhongbo，et al.Study of the Model of Spring Flow Forecast:A Case Study of the Xiaonanhai Spring Catchment［J］.Journal of Sichuan University(Engineering Science Edition)，2011，43(1):77-82.(in Chinese))

［8］张有天.裂隙岩体渗流数学模型研究现况［J］.人民长江，1991，22(3):1-10.(ZHANG You-tian.Current Research in Modeling Seepage Field Fractured Rock Mass［J］.Yangtze River，1991，22(3):1-10.(in Chinese))

［9］王环玲，余宏明.结构面模拟分析在岩体质量研究中的应用［J］.勘察科学技术，2002，(4):25-29.(WANG Huan-ling，YU Hong-ming.Discontinuity Modeling Analysis in the Application of Rock Mass Quality Research［J］.Site Investigation Science and Technology，2002，(4):25-29.(in Chinese))

［10］WHITE W B.Role of Solution Kinetics in the Development of Karst Aquifers［C］∥IAH.Proceedings of the 12th IAH Congress，Huntsville，Alabama，September 21-27，1977:503-517.

［11］袁道先，刘再华，林玉石，等.中国岩溶动力系统［M］.北京:地质出版社，2002.(YUAN Dao-xian，LIU Zaihua，LIN Yu-shi，et al.Karst Dynamic System in China［M］.Beijing:Geological Publishing House，2002.(in Chinese))

［12］GABROVŜEK F，DREYBRODTB W.A Model of the Early Evolution of Karst Aquifers in Limestone in the Dimensions of Length and Depth［J］.Journal of Hydrology，2001，240(3/4):206-224.