典型石漠化地区岩溶水系统循环演化分析*

姚邦杰 刘 琦② 任 标 王 涵

( ①同济大学土木工程学院地下建筑与工程系 上海 200092)

( ②岩土及地下工程教育部重点实验室( 同济大学) 上海 200092)

0 引 言

岩溶石漠化是指在岩溶环境下，受到人类不合理的社会活动和经济活动等的干扰和破坏造成的土壤严重侵蚀、土壤肥力和生产力大幅下降，并导致基岩大面积裸露且出现类似于荒漠化景观的生态环境问题，它是土地荒漠化的主要类型之一(王世杰等，2003)。

岩溶地区因为水对可溶性岩石的侵蚀以及两者之间的能量物质交换形成了该地区特有的地形地貌形态，这种地形地貌形态可以分为地下和地表两个子系统。其中管道、溶洞、地下暗河等组成了地下子系统，而峰丛、溶沟、溶槽、溶隙、岩溶漏斗、天窗、落水洞等是地表子系统的组成部分( Ford，1989; 卢耀如，1999) 。岩溶水径流在地表发生并分配，然后在地下子系统中输送和调蓄。岩溶系统存在孔隙、裂隙和管道3 种介质空间，孔隙流、裂隙流和管道流并存。因此，岩溶地下水在结构上表现为功能上的耗散结构、三维的空间地域结构、二元流场形态结构还有不均一的双重含水介质结构( White，1969) 。岩溶地下水一般不符合达西定律，表现为紊流或者混合流( Ford，1989; 石朋等，2012) ，正是由于岩溶水系统循环演化的特殊性，水不仅会携带土壤发生地表流失，而且会发生地下漏失，导致岩溶石漠化现象突出。岩溶水系统循环演化受到自然和人为因素控制。自然因素包括降雨量、蒸发量、年径流量、地形地貌、地质构造、岩性、水文地质条件等; 人为因素包括土地利用形式、水利水电工程、向河流排放污染水源等( Stringfield et al.，1979; Rosenberry et al.，2008) 。国内外学者进行了大量的研究，其研究方法主要包括:室内实验法、野外调查试验法、理论分析法和数值模拟法( 史运良等，1992; 卢卫中，1995; 袁道先等，1996; 王腊春等，2000; 韩宝平等，2004; 章程等，2004; 刘延惠等，2005; 曾科，2012) 。但目前少有对典型岩溶石漠化地区的岩溶水系统演化问题进行系统性研究，而这一循环过程恰恰是导致石漠化地区土壤地表流失和地下漏失方向和趋势的前提条件和关键性问题。

本文以典型岩溶石漠化地区贵州贞丰—关岭花江高原峡谷区为研究区，对该区的大气降水、河水、泉水进行水化学、氢氧同位素分析，研究岩溶水的补、径、排特征，水化学特征的变化; 利用Modflow对研究区进行地下水数值模拟，用PEST 法( 参数评估) 进行模型识别与验证，通过构建研究区的地下水流场，计算研究区的降雨入渗量、地下暗河的排泄量、地下水泄流量、水均衡情况等，以定量研究区内的岩溶水系统循环演化过程，从而为进一步研究石漠化地区土壤地表流失和地下漏失的方向和趋势提供相关理论依据。

1 研究区概况

研究区位于贵州省安顺市关岭县的北盘江北岸与打邦河西岸所夹的三角地带，东经105°18'45″～105°52'30″、北纬25°33'45″～25°50'00″，研究区内岩溶地貌发育，地形高差极大，地形复杂。研究区年均降雨量约为1100 mm，年均温约18 ℃，海拔376 ～1828 m，研究区河流两岸高差由于河流的下切作用非常大，是典型的岩溶高原峡谷地貌。该区石漠化发育程度强烈，石漠化面积约占研究区总面积的49%，且主要是中度和强度石漠化，已经对人类的生存发展造成威胁。

该区主要出露三叠系、白垩系及第四系地层，其中三叠系地层出露面积占95%以上，碳酸盐岩出露面积约占85%，根据岩性、岩层的水理性质及地下水赋存状态，将研究区地下水划分为3 类，包括碳酸盐岩裂隙溶洞水、碎屑岩碳酸盐岩溶洞裂隙水和层状岩类裂隙水。研究区内共有主要地下河两条，岩溶管道总长度约39 km，岩溶大泉20 个。研究区主要发育有4 组裂隙，并且主要发育1 组北西向的裂隙，岩溶发育的主控方向即主要受北西向裂隙控制，几条主要地下暗河走向均为北西向，主要分布于研究区打邦河西岸，各含水类型及地下河的分布见图1。

2 岩溶水化学特征及补、径、排特征

2018 年9 月，在研究区进行野外调查并采集雨水样4 组，在采集雨水样2 d 后采集泉水样20 组，河水样2 组( 采样点见图1) ，基本覆盖研究区泉水出露点，采样后进行氢氧同位素和水化学成分测定，结果见表1。

图1 研究区岩溶水系统水文地质图Fig. 1 Hydrogeological map of karst groundwater system in the study area

因随着海拔的增加，气温相应降低，此时由于蒸发分馏作用的影响，大气降水中氢氧同位素组成会发生线型变化，因此研究区大气降水中氢氧同位素的组成也具有显著的高程效应( 陈中笑等，2010) ，据此可推算出地下水的补给高程，它是判断地下水的补给位置以及补给范围的基础，岩溶地区地下水循环更替较快，该方法更加有效( 黄荷等，2016) 。通过研究区水样分析计算结果表明，花江雨水中δD值与海拔成正相关关系: H=－28.976 δD－998.25( 图2) 。每个主要岩溶泉的氢氧同位素值代入研究区雨水δD 值与降雨高程的关系公式，计算出主要岩溶泉的大致补给高程，并根据结果标注到地图上( 图3) 。大部分泉水的补给范围与岩溶发育主控方向( 北西—南东) 一致，泉水的补给区域主要受断层、地下管道控制，其中16 号泉水( 地下暗河) 补给范围最大，从西北边的上关镇到东北部的坝坎周围均为其汇水范围。根据泉水补给范围推测出研究区地下水流向主要为北西—南东向，靠近南部北盘江时地下水流向转为近北—南向，局部地区如花江镇东南部地下水流向因地形高差原因会有一定变化。

3 研究区地下水数值模拟

为了深入了解研究区岩溶水资源转化量的关系，本文使用MODFLOW 软件进行地下水数值模拟，根据室内外资料构建地下水流场，计算研究区的水均衡情况、降雨入渗量、地下河管道的排泄量、地下水泄流量等，以定量研究区内岩溶水转化关系。

表1 研究区泉水、地表水、雨水同位素及水化学分析数据Table 1 Spring，surface water，rainwater isotope and water chemical analysis data of the study area

图2 研究区雨水δD 值的高程效应Fig. 2 Elevation effect of rainfall δD value in the study area

3.1 地下水数值模拟

按照研究区的水文地质条件圈定研究区的边界，并将研究区的所有汇水面积圈定为模型范围，最终模型横向长31.7 km，纵向长38.3 km。其中大面积发育的碳酸盐岩为模型的主要含水介质。研究区的边界主要分为河流和断层两类，北部由两条大型压性断层构成边界，南部由北盘江和打邦河两条大型河流构成研究区边界。模型的边界条件主要分为4 类，第1 类为大型河流构成的河流边界，分为两段，模型南面为北盘江，东面为北盘江支流打邦河，大部分地下水补给入这两条河流，因此将这两条河流概化为模型的河流边界和流出边界; 第2 类为地下暗河构成的排水沟边界，模型边界内主要发育有许凹地下河，运用MODFLOW 中的“排水沟”子程序来模拟该条地下河，为流出边界; 第3 类为模型西南边海拔最高处的定流量边界; 第4 类为补给边界。在稳定流计算时该地区的降雨补给量设定为400 mm·a－1。非稳定流计算时，降雨补给量按气象站实测值的0.33 倍给出。模拟范围、研究区边界及模型边界条件如图4 所示。

图3 研究区主要岩溶泉推测补给范围Fig. 3 Speculation range of major karst springs in the study area

表2 水文地质参数初值Table 2 Initial values of hydrogeological parameters

依据上述水文地质概念模型，可将研究区的裂隙岩溶地下水流问题概化为非均质各向异性介质中的三维非稳定流问题，研究区在垂向上分为6 层，含水层均为潜水含水层，水文地质参数初值如表2 所示。对模拟区进行网格剖分。模型在垂向上除了原有的分区以外不进行更细致的划分。模型的网格剖分总体上采用154×159 m2的形式。剖分后活动单元格数约为24 500 个，模拟面积1214.1 km2。以研究区内主要岩溶泉的高程值作为该点地下水位的拟合值，采用MODFLOW 中的PEST 算法进行参数反演。

3.2 数值模拟结果

3.2.1 地下水位变化

研究区地下水流场模拟结果如图5 所示，地下水流向大致呈北西—南东向，分别汇入南部北盘江和东部打邦河，地下水流向与该地区岩溶发育主控方向一致，较好反映了研究区的地下水流场。OW-4观测点的水位变化情况( 图6) ，该区地下水位年变化幅度较大，季节性强，3 ～5 月为地下水枯水期，8～10 月为地下水丰水期，地下水位变化情况相对降雨量变化情况有一定的滞后性。

3.2.2 地下暗河流量

图4 模型范围、研究区边界及模型边界条件Fig. 4 Model range，study area boundaries，and model boundary conditions

模拟得到地下暗河流量变化情况( 图7) ，地下暗河流量变化与降雨量变化情况存在一定的滞后性，地下暗河流量随雨季来临而上升速度加快，模拟地下暗河流量的回落速度在降雨量减小时较慢。结合模拟地下河出口的总排泄流量与实际观测流量的拟合效果分析，模拟雨季最大流量小于实测雨季最大流量，而模拟枯季最小流量大于实测枯季最小流量，可见模拟流量变化情况比实际情况稍缓，但误差不大，模拟过程和实际条件下地下水系统排泄量的变化特征基本吻合。

3.2.3 水均衡

模拟计算出的研究区年均地下水均衡情况( 表3) ，本次模拟期中降雨补给量是主要的补给来源; 侧向补给量是次要补给来源。河流侧向排泄量占排泄量的60.21%，虽然占比较大，但唯一一条地下河的管道排泄量占总排泄量的24.92%，符合岩溶地区的实际情况。

将2018 年每个月研究区的水均衡情况归纳如图8 所示。在雨季，大量增加的地下水无法迅速通过地下暗河或河流排泄，其中的大部分被地下含水层储存，而地下水向地下暗河及河流的排泄量迅速增大。在旱季，因该地区几乎没有大气降水的补给，地下含水层中储存的地下水开始发挥补给功能，而地下水向地下暗河及河流的排泄量迅速减小。一年之中，侧向排泄量和侧向补给量变化较小，地下暗河的排泄量变化幅度显著大于河流排泄量的变化幅度。

根据模拟计算结果看出，研究区33%的大气降水通过表层岩溶裂隙带的下渗作用转化为地下水，其中24.9%的地下水通过地下暗河排泄、60.2%的地下水通过河流排泄作用进入北盘江和打邦河成为地表水，其余的14.9%的地下水侧向排泄成为区域外的地下水。

4 结 论

通过野外调查、水样分析及数值模拟等方法，对贵州贞丰—关岭花江高原峡谷石漠化治理示范区岩溶水系统的补、径、排特征，水化学特征及水资源量等进行分析，得到如下结论:

图5 OW-4 观测点地下水位模拟值Fig. 5 Groundwater level simulation value of OW-4 observation point

图6 模拟地下水流场和地下水位拟合点Fig. 6 Simulated groundwater flow field and groundwater level fitting point

表3 年均地下水均衡表Table 3 Annual average groundwater balance table

(1) 研究区内大气降水的氢氧同位素组成具有显著的高程效应，据此计算出了研究区主要岩溶泉的补给高程。综合水化学分析推断出的补给区地层，以及氢氧同位素分析计算出的补给高程，得到了研究区主要岩溶泉的补给范围。得出研究区河水上下游的水化学、氢氧同位素变化与河流两岸的海拔、岩性的变化高度一致，可见地下水与河水之间联系密切。

图7 地下暗河流量变化模拟情况Fig. 7 Simulation of underground dark river flow changes

图8 研究区2018 年水均衡变化情况Fig. 8 Water balance changes in the study area in 2018

( 2) 建立了关岭—贞丰花江地区MODFLOW 数值模型，根据实际情况将含水岩层分为6 层，模拟过程和实际条件下地下水系统排泄量的变化特征基本吻合，根据模拟结果得到，研究区33%的大气降水通过表层岩溶裂隙带的下渗作用转化为地下水，其中24.9%的地下水通过地下暗河排泄、60.2%的地下水通过河流排泄作用进入北盘江和打邦河成为地表水，其余的14.9%的地下水侧向排泄成为区域外的地下水。其结果符合岩溶地区的实际情况，但由于岩溶地区地下水流情况十分复杂，且地层资料精度较差、地下水位观测数据缺乏，地下水的模拟仍有一定误差。

分析结果基本反映了实际岩溶水系统排泄量变化及地下水流向特征，初步阐明了该区域岩溶水系统循环演化特征，为岩溶石漠化地区的土壤地表流失和地下漏失的方向和趋势研究及防治提供理论依据。