基于MODFLOW-CFP 的贵州大井流域岩溶地下水数值模拟

基于MODFLOW-CFP 的贵州大井流域岩溶地下水数值模拟

党志文  1 ，邵景力  2 ，崔亚莉  2 ，李 军  1 ，宫志强  3 ，赵良杰  4 ，梁永升 5

（1.河北建筑工程学院, 河北省水质工程与水资源综合利用重点实验室, 河北 张家口 075000；2.中国地质大学水资源与环境学院, 北京 430074；3.河北省地质资源环境监测与保护重点实验室, 河北省地质环境监测院,河北 石家庄 050022；4.中国地质科学院岩溶地质研究所/自然资源部、广西岩溶动力学重点实验室,广西 桂林 541004；5.河北省地矿局第三地质大队, 河北 张家口 075000）

摘 要： 贵州大井流域岩溶分布广泛，岩溶水是当地人民生产和生活的主要来源。由于对岩溶水资源的不合理开发利用，水资源短缺现象经常发生。大井流域水文地质条件复杂，管道−多孔介质双重介质特征明显。文章采用MODFLOW-CFP 耦合模型对大井流域展开数值模拟，进而掌握大井流域地下水运动规律、准确评价岩溶水资源，促进其合理开发利用。结果表明：大井流域管道与多孔介质交换量为6 719.1 m 3 ·a −1 ，主要集中在上游和中游；总补给水量为10 977.3×10 4  m 3 ·a −1 ，补给模数为133.495 m 3 ·km −2 ·a −1 ，其中降雨汇入量和降雨入渗量占总补给量的81.35%，而总排泄量为10 813.47×10 4  m 3 ·a −1 ，主要在地下河出口排泄。

关键词： 大井流域岩溶水；MODFLOW-CFP 耦合模型；岩溶管道模拟

0 引 言

岩溶地区含水介质具有多重复杂性，数值模拟难度较大，故其发展十分缓慢。传统的黑箱模型和概念模型物理意义十分局限，无法反映含水系统内部完整的水文过程 [1] 。分布式模型近几年得到长足发展，它基于水文动力学机制，能较好模拟岩溶水系统内外的水文过程。目前国内外比较流行的有等效多孔介质模型、管道模型和多孔介质−管道模型 [2−5] 。

等效多孔介质模型将岩溶含水系统概化为均质模型，并且假定地下水的流动呈层流形态，符合达西定律。陈皓锐等 [6]  利用多孔介质模型构建沧州市吴桥县的三维地下水流动模型。Dong Y H 等 [7] 对多孔介质模型进行改善，提高了模拟的精度。等效多孔介质模型不需要管道系统的详细数据，模型建立相对容易，但其概化的方法导致无法反映出不同岩溶含水介质差异的特点，且无法模拟不同介质之间的水力交换情况 [8] 。管道模型忽略除管道外的其他介质，将岩溶含水系统作为独立的管道网络，地下水在管道中呈现紊流或者层流状态，取决于雷诺系数，当管道流为层流时，采用Hagen-Poiseuille 公式计算，当管道流为紊流时，采用Darcy-Weisbach 公式计算 [9] 。Springer 将管道模型应用于以管道流为主的岩溶含水系统。SWMM 模型由美国环保署研发，采用圣维南方程模拟管道流，该模型后来被吴月霞等 [10] 相继引入中国管道系统的水文模拟中。王益伟、赵良杰、薛亮、梅向阳等 [11−14] 对管道系统的几何参数进行了探讨，并对敏感性进行了分析，张春艳等 [15] 利用自行设计的物理模型对落水洞水位与水位情景响应做了深入研究。杨郑秋、钱家忠、刘丽红、李向全等 [16−19] 在MODFLOW 源代码增加管道紊流程序，将管道流模块CFP 耦合到MODFLOW 模型中，提出了MODFLOW-CFP 耦合模型 ，该模型基于模块有限差分，在理论上较为成熟。管道模型能真实反映岩溶含水系统的管道流特点，但模型建立需依托大量管道数据，在一些勘探难度较大的地区难以应用。多孔介质−管道模型根据介质将岩溶含水系统分为两部分 [20] ，一部分为多孔介质、孔隙组成的多孔介质系统，另一部分为管道介质系统，两部分独立运行，在二者之间存在水流交换机制。

贵州大井流域属于典型岩溶峰丛洼地，洪涝灾害频发，在丰水季节，流域的地下河下游时常遭受大规模淹没，而在枯水季节地下河上游出现不同程度的缺水 [21] 。洪涝灾害的频发和水资源无法得到合理利用严重制约着当地经济的发展。另外，当地人的生产生活用水主要依赖表层岩溶泉，导致表层岩溶泉的流量逐年减少。为解决这些问题，掌握贵州大井流域的地下水运动规律十分必要。本文针对研究区岩溶含水介质二维特征明显的特点，采用MODFLOW-CFP 耦合模型进行数值模拟，建立大井流域地下水流模型，探讨研究区的地下水运移机制，掌握两种含水介质之间的水量交换规律，以期为未来地下水资源的可持续利用提供依据。

1 研究区概况

1.1 大井流域概况

研究区位于贵州高原南部斜坡地带，属惠水县、平塘县、罗甸县管辖，面积为82.23 km 2 。地势总体特征表现为北高南低，地形起伏强烈，相对高差250～400 m，最大相对高差达600 m。东部及东北部海拔相对较高，在1 200 m 以上，最低点位于地下河出口位置，高程为430 m。大井岩溶流域地处亚热带湿润季风气候区，干湿季节分明，总体气候特征表现为：降水量、相对湿度由北向南逐渐递减，气温从北向南逐渐增高，多年平均气温为19.6 ℃，多年平均蒸发量为1 253.8 mm，多年平均降水量为1 178.2 mm。每年5-9 月为集中降水期，降水量占全年降水总量的76%。

研究区地处黔中山原向广西峰林平原过渡的斜坡地带，岩溶地貌以峰丛洼地为主，出露地层从二叠系至三叠系均有分布。地层岩性以碳酸岩为主，三叠系的凉水井组和小米塘组是研究区的主要含水层，在整个研究区几乎都有发育；研究区的西部及东南部有三叠系的大冶组分布，成分以灰岩为主，但由于掺杂着泥岩，表现为弱透水性；研究区的东北边界发育二叠系的吴家坪组石灰岩，是主要隔水层。其在构造上位于扬子准地台黔南台陷之贵定南北向构造变形区，特有的构造应力场造就了区域构造框架中分布最广、规模最大的经向构造体系 [21] 。大、小井出口是经向构造体系和纬向构造体系的复合部位，研究区从地下河入口到出口发育一条南北向的董当张性断裂，在两条地下河交汇处发育NE-SW 向的节理，这两条大的构造控制岩溶发育方向，同时在研究区碳酸盐岩分布的地方岩溶也较为发育（图1）。

图1 大井流域水文地质概图 Fig.1 Hydrogeological survey of Dajing basin

1.2 大、小井间的水力联系

除了研究区内部2 条地下河有密切水力联系外，左侧地下河还与小井流域的地下河有复杂水力联系，主要体现为以下两点：从研究区内航龙入口处大井地下河与小井地下河开始分开，但在分离点处大井与小井的水量分配一直没有弄清楚；由于地下河形态多样、暗河较多，大小井在边界上仍存在一定的水力联系。

针对第一点问题，经过实际调查发现：在丰水季节，航龙入口西侧存在一条季节性河流，约30%的入口水量通过这条河流流入小井流域，剩余70%全部转为暗流流入大井地下河；在枯水季节，入口流量全部转化为暗流流入大井地下河。针对上述第二点问题，贵州省地质调查院在2003 年开展了大、小井地下河系统示踪试验，结果表明：在航龙至马鞍寨区间，大、小井地下河水力联系密切，而在马鞍寨之后两条地下河几乎没有水力联系，相互独立（图2）。

图2 两条地下河间的水力联系(来源：贵州省地调院，2003) Fig.2 Hydraulic connection of two underground rivers(Source:Guizhou Provincial Geological Research Institute, 2003)

2 水文地质概念模型

2.1 含水层结构

结合大井流域的水文地质情况，将研究区从地表到基岩底板概化为潜水含水层。含水层的顶板为地面标高（通过DEM 数据获得），底部边界为基岩底板。

2.2 边界条件

大井流域是一个相对完整和独立的水文地质单元，大井流域的具体划分是以地下水、地表水的分水岭和研究区的构造（断裂、节理等）对地下水的补给、径流和排泄条件所起的相对控制作用为基础。研究区的西部是大井流域与小井流域的分水岭，区域的东部和东南部边界为山区，海拔较高，形成分水岭。故将东部和西部边界概化为隔水边界。研究区东北部沿着航龙、旧宅、塘边一带为二叠系吴家坪组，岩性为渗透性很差的泥岩和页岩，处理为隔水边界。北部航龙附近是伏流入口，处理为流入边界。南部为大井出口，处理为流出边界。

2.3 水文地质参数

2.3.1 多孔介质

研究区多孔介质含水层参数依据水文地质资料和前人工作成果，并在模型识别过程中进行修正，得到降雨入渗补给系数、渗透系数、给水度（表1）。蒸发较弱，平均蒸发系数0.1，蒸发极限埋深与第四系厚度一致，取平均值2 m。

表1 校正后模型水文地质参数 Table 1 Corrected model hydrogeological parameters

2.3.2 管道

在CFP 模块中需要提供的管道参数包括管道尺寸、管道曲率、管道粗糙系数、管道壁与多孔介质的交换系数、层流紊流互相转化的高低雷诺数以及地下水水温。根据前人报告和钻探数据获得地下河管道平均直径、摩擦系数和上下界雷诺数、管道曲折程度等。结合示踪技术和模型运行结果反演管道的水文地质参数，得到管道平均坡度为2.5%，根据地下河入口和排泄出口的高程，以及落水洞、洼地的高程计算出管道中每个节点高程，管道曲折度的数学定义为长度与管道起始点直线距离的比值，利用ArcGIS 功能中的距离测量可间接计算出管道曲度，其平均曲度为1.05。

2.4 地下水补给与排泄

研究区地下水补给项包括大气降雨补给、伏流入口流入，研究区北部的侧向流入补给。大气降雨补给包括降雨入渗和降雨过后的汇流流入。排泄项主要包括地下河出口排泄、蒸发和泉排泄。

蒸发排泄和泉排泄通过前期的调查和监测获得；大气降雨数据在中国气象站下载获得，降雨量与降雨入渗补给系数的乘积即为降雨入渗量，再减去地表蒸发量即为汇流流入量；同时在本次研究之前，已完成大、小井流域SWAT 地表水模型的建立，在该模型基础上，获取研究区北侧的流入补给量以及汇水总量在各个小流域上的分配。

3 地下水流数值模型

依据地下水动力学原理，在遵循质量守恒和能量守恒的基础上，利用软件对大井流域地下水进行数值模拟，通过模型识别和水均衡分析，验证模型的准确性。针对研究区的含水介质表现为管道−多孔双重介质的特点，本次数值使用MODFLOW-CFP 软件搭建多孔介质−管道耦合模型。

3.1 模型输入

模型的输入包括含水层结构、边界条件以及水文地质参数，在前文水文地质概念模型中已介绍，只需将其编辑输入到MODFLOW-CFP 耦合模型中即可。

3.2 时空离散

该模型在垂直方向分为一层，水平方向上结合研究区面积大小和计算复杂程度，将研究区剖分成180 行×120 列的矩形网格，每个网格代表实际大小100 m×100 m 的区域。在管道位置，采用CFP 模式中的CFPM2 模块刻画，节点仍采用100 m×100 m 的均匀剖分格式，共计299 个节点，301 条管道流。

选取2014 年为模拟期，以月为应力期进行模拟。

3.3 模型识别和验证

3.3.1 流场形态

由于研究区水文地质资料不足，无法直接获得初始流场。根据已有DEM 数据、降雨资料、大井流域流入和流出数据，以及研究区各套岩性的水文地质参数的经验数据，模拟2014 年的稳定流，获得稳定流初始流场（图3）。由于采用的降雨数据、流量数据等是月平均值，该稳定流场代表的近似于平水季节的流场形态。流场主要受地形地貌和地下河管道控制，远离地下河地区的流场与地势起伏近似，靠近地下河的地区管道特征明显，即管道位置流场背离地下水流方向。

图3 初始流场 Fig.3 Piper diagram of initial flow

图4 是模拟区域枯水季流场形态，地下水位整体上从北向南逐渐降低，基本上与地势起伏一致，东北部地下水位较高，水位最高值接近880 m，在大井地下河出口附近水位最低，最低值为420 m。管道附近水位略低于周围的水位，周围多孔介质地下水流向趋向管道水流方向，在中游附近尤为明显。研究区自上游至下游等值线由稀疏逐步变密集，说明中上游地势起伏小、下游地势起伏大。

图4 枯水季流场图 Fig.4 Piper diagram of flow in dry season

相对于枯水季节，丰水季整体水位明显抬升。如图5 所示，最高水位为883 m，同样位于东北地势较高处，最低水位为424 m，在地下河出口，但明显高出基准面4 m 多，说明在丰水季节，出口附近被大面积淹没。另外，管道特征更为明显，管道附近水位明显低于周围水位，在中游附近管道水位与两侧水位相差接近10 m，体现出中游附近地势切割强烈且丰水季陡坡急流的特点。

图5 丰水季流场图 Fig.5 Piper diagram of flow in rainy season

3.3.2 水位拟合

研究区水位监测数据有限，且缺少一个完整年的水位监测，故在研究区地下河入口处选取丰水季和枯水季的一段时间进行拟合。图6 是2014 年1 月和7 月的水位拟合情况，从图上可知，1 月份观测孔水位动态曲线整体拟合效果较好，观测水位与模拟水位在数值上基本一致，且能反映出地下水的动态变化规律；7 月份的水位模拟数据与观测数据相比，地下水动态趋势拟合较好，对降雨的动态反应略为滞后，但观测水位与模拟水位数值上有点偏差，模拟水位整体高出观测水位2～3 m。分析其原因：可能是由网格大小和观测孔附近地势起伏大引起的，单元格大小是100 m×100 m，而在该范围内与观测孔所在位置的最大相对高差达到10 m 以上，导致单元格代表的平均数值与实际观测数值存在出入。

图6 水位拟合图 Fig.6 Piper diagram of water level fitting

3.3.3 出口流量拟合

图7 是大井流域2014 年1-12 月的出口流量拟合情况，从过程线的总体拟合情况来看，观测流量和模拟流量在数值和变化趋势上基本一致，基本反映出大井流域下管道排泄量月变化特征，4-9 月份出口流量较大反映出丰水季节的特点，10-3 月份出口流量明显减小反映出枯水季的特点。模拟流量与实测流量在趋势上基本保持一致，在丰水季节，模拟值略小于实测值，枯水季节模拟值略大于观测值。

图7 出口流量拟合图 Fig.7 Piper diagram of outlet flow fitting

3.4 模拟结果

3.4.1 滞留时间

图8 是CFP 模型模拟出的地下河管道从上游到下游滞留时间的分布。滞留时间在上游的前11 处管道急剧减小，前11 处管道与多孔介质的交换水量较大，导致短时间内通过管道的流量迅速增大，从而滞留时间急剧减小。在第63、101 和171 管道处也会看到滞留时间迅速下降的趋势，反映的是落水洞位置处的降雨集中补给。进入中游，滞留时间缓慢下降，反映多孔介质一直对管道的补给，但补给量相对上游较小。进入下游特别是接近大井出口的位置，滞留时间不再发生变化，管道与多孔介质的水位几乎持平，二者不存在水量交换。

图8 管道滞留时间图 Fig.8 Piper diagram of pipe’s residence time

3.4.2 平水期管道与多孔介质交换量

为进一步研究从上游到下游管道与多孔介质的交换关系，选取代表平水期的部分管道进行研究，如图9 所示，交换量代表多孔介质对管道的日补给量，交换量整体上呈现逐渐减少的趋势，在前109处管道，交换量是正值，代表多孔介质对管道的持续补给；109 处管道开始到114 处管道出现了短暂的负交换量，表示在这一小段地下河出现管道补给多孔介质的现象，这可能是由于该处管道直径突然较小，引起管道内压强突增，进而出现逆补给现象；从114 处管道以后，交换量几乎为0，表示管道与多孔介质之间不再有水力交换发生，也表征两种介质的水位开始保持一致。从以上分析更进一步说明在地下河的上游、中游管道中的水未处于饱和状态，不断得到补给，到了下游管道开始充满水，不再接受多孔介质补给。

图9 从上游到下游管道与多孔介质交换量变化图 Fig.9 Diagram of changes in pipeline and porous media exchange capacity from upstream to downstream

3.4.3 水均衡

多孔介质−管道模型整体的均衡情况见表2，本次模拟期中降雨补给量为8 930×10 4  m 3 ·a −1 ，是模型主要补给量，占总补给项的81.36%，其中，管道降雨入渗，即降雨通过落水洞进入管道的水量为3 968.9×10 4  m 3 ·a −1 ，多孔介质降雨入渗量为4 961.2×10 4  m 3 ·a −1 。北 部 航 龙 入 口 的 补 给 量为2 047.2×10 4  m 3 ·a −1 ，占18.65%。排泄项中地下河排泄量为10 688×10 4  m 3 ·a −1 ，占总排泄项的98.84%，是最主要的排泄项，与实际情况中地下河流速大的特点相符，同时这个排泄量几乎等于降雨总量。排水沟指的是东北部狭长的低洼地带排入到地下河中的水量，为114.47×10 4  m 3 ·a −1 ，占总排泄项的1.06%。蒸发量相对较少，这与大井流域地下水的埋藏较深有关，仅在地下河入口、出口以及低洼地带存在部分蒸发。另外研究区还存在29 处下降泉，排泄量较少，仅为6.5×10 4  m 3 ·a −1 。

表2 地下水均衡分析表 Table 2 Groundwater balance analysis

整个模型的均衡情况能体现出整体的源汇项，但无法体现管道与多孔介质之间的交换情况。而无论在实际情况中还是在模拟过程中，交换量十分重要，需要单独计算管道与多孔介质的均衡情况。

单独计算管道模型的均衡情况如表3，管道中补给项主要由降雨和管道交换组成，其中降雨项指大气降雨通过落水洞的直接补给，水量为3 968.9×10 4  m 3 ·a −1 ，占总补给项的37.13%，管道交换项即多孔介质对管道的补给，水量为6 719.1×10 4  m 3 ·a −1 ，占总补给项的62.87%。排泄项全部为管道排泄，排泄总量达到10 688×10 4  m 3 ·a −1 。

表3 管道水均衡分析 Table 3 Pipe’s groundwater balance analysis

单独计算多孔介质的均衡情况如表4，多孔介质中补给项主要由降雨和北部河流流入组成，其中降雨项指大气降雨对多孔介质的面状补给，水量为4 961.2×10 4  m 3 ·a −1 ，占总补给项的70.79%，北部航龙入口流入水量为2 047.2×10 4  m 3 ·a −1 ，占总补给项的29.21%。排泄项由管道交换、地下河流出、蒸发以及泉排泄组成。其中多孔介质对管道的补给量为6 719.1×10 4  m 3 ·a −1 ，占总排泄项的98.17%，排水沟项指的是研究区东北部低洼地带排泄到两条地下河的水量，水量为114.47×10 4  m 3 ·a −1 ，占总排泄项的1.67%，蒸发和泉排泄量较小，占总排泄项的0.16%。

表4 多孔介质水均衡分析表 Table 4 Fracture water balance analysis

4 结论与建议

以贵州大井流域为研究区，通过MODFLOWCFP 建立管道、多孔介质双重介质模型，结合流场形态、观测孔水位拟合定性和定量地验证模型的准确性，得出以下结论：

（1）研究区总补给水量为10 977.3×10 4 m 3 ·a −1 ，补给模数为133.495 m 3 ·km −2 ·a −1 ，其中降雨汇入管道量为3 968.9×10 4 m 3 ·a −1 ，占总补给量的36.16%，降雨入渗量为4 961.2×10 4 m 3 ·a −1 ，占总补给量的45.20%，地下河入口流入量为2 047.2×10 4 m 3 ·a −1 ，占总补给量的18.65%。总排泄量为10 813.47×10 4 m 3 ·a −1 ，其中地下河 出 口 流 出 量为10 688×10 4 m 3 ·a −1 ，占 总 排 泄 量的98.84%，泉和蒸发排泄量占了剩下的1.16%。从大井流域的水量均衡情况看出，降雨量为主要的补给来源，地下河出口流出量为主要的排泄去处，二者几乎相等；从大井流域与外界的联系来看，外界对大井流域的补给量即大井流域流入量远远小于大井流域流出量，这反映出大井流域的水资源利用程度偏低；

（2）研究区管道与多孔介质介质的交换量为6 719.1 ×10 4 m 3 ·a −1 ，几乎全部是多孔介质补给管道。另外，无论从滞留时间分析，还是从地下河上游到下游的两种介质间的交换量分布，都表明地下河上游管道与多孔介质水量交换频繁，下游几乎没有水量交换，即从上游到下游交换量逐渐下降，而在某些管段位置会发生下降幅度增大的现象，反映出该处管段处在洼地或落水洞，降雨过后的汇水使管道内水位充满，多孔介质对管道的补给量进而减少；

（3）大井流域虽然降水充沛，但水资源短缺现象仍时有发生，主要与其开发利用方式有关。根据上述分析，提出两点开发利用方式取代传统的过度利用岩溶泉的取水方式：一是增加地下河出口的引水，二是截取上游、中游的汇水量作为备用。