考虑山坡和洼地水力联系的半分布式喀斯特水文模型——以陈旗流域为例

薛冰贤，陈 喜,2，陈 曦，张志才，程勤波

(1. 河海大学 水文水资源与水利工程科学国家重点实验室, 南京 210098；2. 天津大学 表层地球系统科学研究院, 天津 300072)

我国南方喀斯特地区有着独特的地貌和水文特征，山坡土层薄、持水能力弱，植被以林地、灌丛为主；洼地土层较厚，多为农田，因此，山坡和洼地滞留水量是植被耗水、农业灌溉用水的重要水源。山坡、洼地土层以下溶沟、溶隙和管道发育，是水土流失的重要通道；暴雨时山坡水流向洼地汇集，极易形成陡涨陡落雨洪过程，引发洪涝灾害。因此，喀斯特地区坡、洼系统水文过程的模拟与预测对生态保护、水资源利用和洪涝灾害防治具有重要意义。

喀斯特地区不同大小裂隙和管道的导水介质导水性和水流特征显著不同，因此地下水流模拟通常采用双重介质模型，如在孔隙介质地下水渗流(基质流)模型基础上增加管道流计算，常用的有MODFLOW-CFP[1]模型；另一类是把导水裂隙概化为等效多孔介质[2]，并考虑导水介质各向异性模拟喀斯特流域地下水动态或水文循环过程，这类分布式模型通常应用于具有详细岩溶裂隙分布、结构、导水性等观测资料的实验流域。集总式、半分布式水文模型通常在汇流计算层次考虑不同介质对水流的调蓄作用[3]概化水流特征，如把蓄水体概化为两个或多个水箱[4]，建立基于落水洞的岩溶半分布式水文模型[5]，模拟和预测降雨-径流响应特征，但对不同介质产流入渗机制和不同地貌单元(如山坡与洼地)以及它们之间水力联系的模拟功能还不足，且喀斯特地区入渗补给与非喀斯特不同，具有表层岩溶带裂隙以及漏水洞、漏斗等快速入渗和水流运动通道。因此，发展兼具资料条件和喀斯特地貌特征的半分布式模型具有重要意义。

本文考虑峰丛喀斯特地貌特征，将流域划分为山坡、洼地两个单元，并考虑裂隙和管道的产流和入渗机制以及饱和带快速流和慢速流的交互作用，构建半分布式水文模型。利用流域出口断面流量以及洼地地下水位等观测资料进行模型参数率定和验证，分析坡、洼系统多重介质径流量分布，为峰丛喀斯特地区水文过程模拟和水资源量评估提供支撑。

1 模型与方法

将峰丛喀斯特流域划分为山坡(H)、洼地(D)两个单元(面积分别为AH和AD)，每个单元内根据蓄水、导水能力不同，划分为具有基质流特征的土壤与细小裂隙区(所占面积比例为α)以及具有管道流特征的大裂隙管道(漏斗)区(所占面积比例为(1-α))；垂向划分为非饱和带和饱和带，进行单元内产流、入渗补给和蓄泄关系以及两个单元之间水力联系计算(图1)。

1.1 非饱和带产流和入渗补给计算

在具有基质流特征的土壤与细小裂隙区，与南方非喀斯特地区产流机制相似，采用新安江模型蓄水容量曲线[6]计算时段产流R。由于喀斯特裂隙发育深，R一部分向下入渗补给饱和带，即IS=ks×R(ks为慢速水箱补给系数)；剩余部分(R-Is)侧向汇至周边大裂隙和管道(漏斗)区。

在具有管道流特征的大裂隙管道(漏斗)区，降落在该区的降雨量P扣除少部分损失后(aP，a为折扣系数)形成的自由水(产流)，加上上述土壤与细小裂隙区部分汇入的径流(R-Is)进入管道；并进行管道蓄量VF,t调蓄计算，管道中一部分蓄量向下补给饱和含水层，即IF=kf×VF,t(kf为快速水箱补给系数)。当管道蓄量VF,t大于其最大蓄量VM时，管道溢出产生地表径流Rs，即Rs=VF,t-VM；否则Rs=0。

1.2 饱和带蓄泄关系演算

基质区和管道区入渗补给饱和含水层中慢速、快速水箱，进行蓄泄演算。在两个单元之间，山坡(H)慢速、快速水箱出流分别流向洼地(D)慢速、快速水箱，再经过洼地调蓄流入地表或地下河出口。模型结构如图1所示。

山坡(H)、洼地(D)单元慢速或快速水箱蓄泄演算如下。

对于山坡(H)单元，水箱j(j等于F或S，代表快速或慢速水箱)蓄量VHj计算公式为：

(1)

对于洼地(D)单元，水箱j蓄量VDj计算公式为：

(2)

任一山坡或洼地单元i，EXi计算公式[7]为：

(3)

式中：ke为水箱交换系数。

假定任一单元i(i等于H或D，代表山坡或洼地单元)中水箱j的出流Rij与蓄量Vij成线性关系，即：

Rij=ηij×Vij/Ai

(4)

式中：ηij为单元i中水箱j出流系数。

Eg采用阿维里扬诺夫公式[8]计算，并考虑地下水埋深空间分布的非均匀性，通常假定服从伽玛分布[9]，Eg计算公式为：

(5)

式中：m为计算单元分区数；Dm为最大平均地下水埋深，m；D0为初始平均地下水埋深，m；Dmax为潜水蒸发极限埋深，m；Γ(γ)为伽玛函数；γ、β分别为伽玛函数的形状、尺度参数；kc为蒸发折算系数；Ep为流域潜在蒸发量，mm。

洼地地下水平均埋深D计算公式为：

D=(VDS+VDF)/AD/sy

(6)

式中：sy为给水度，本文取0.05。

模型应用于西南典型喀斯特小流域，面积很小，因此忽略径流的汇流过程，洼地地表径流Rs以及慢速、快速水箱出流Rij进入河道从流域出口断面流出。

2 流域概况及资料

陈旗流域位于黔中丘原盆地区的普定岩溶盆地(26°15′36″ ～26°15′56″、东经105°43′30″ ～105°44′42″)，属于亚热带季风湿润气候区，流域面积0.9 km2。利用DEM数据将流域划分为山坡、洼地两个单元，面积分别为0.73 km2、0.17 km2。流域多年平均降水量1 300 mm，降雨主要集中在5-8月份，年均气温15.1 ℃。具有贵州典型的峰丛山体、峰丛洼地地貌及喀斯特水文地质特征[10]。

流域内在阳坡、火烧坡位置(图2)分别设立了HOBO小型气象站，获取实时降雨、气温、气压、风速、湿度、太阳辐射等气象资料，利用彭曼公式计算潜在蒸散发量。流域出口处修建了三角堰并分别放置了HOBO水位自计仪观测水位，获取5 min间隔的实时水位数据，通过堰流公式计算获得流域的流量数据。利用流域洼地内4眼观测井的水位数据，推求流域洼地地下水埋深。

图2 陈旗流域图Fig.2 Chenqi watershed

模型率定期为2016年10月8日至2017年10月30日，历时9 299 h；验证期为2017年10月30日至2018年6月12日，历时5 406 h；计算步长1 h。

3 参数率定

3.1 目标函数

目标函数是为了确定模拟过程与实测过程的吻合程度，首先分别采用流域出口流量和洼地地下水埋深的纳什效率系数NSE[11]为单目标函数，再结合流量和地下水埋深的纳什效率系数构建多目标函数[12]，计算公式为：

(7)

式中：σo,Q、σo,D分别为实测流量Q、地下水埋深D的标准差；n为数据时序长度；NSEQ、NSED分别为Q、D的纳什效率系数。

3.2 参数敏感性分析与优化

采用全局敏感性分析方法，分析各个参数以及它们之间相互作用对模型输出的影响。采用Monte Carlo法随机生成各个参数大量的样本[13](本文生成50000组样本)，利用蒙特卡洛分析工具箱(MACT)中的RSA方法对参数和目标函数进行归一化处理，并计算和绘制累计频率分布图[14]，利用频率分布图计算参数敏感性指数(见表1)。敏感性指数大于0.1时，参数敏感；敏感性指数小于0.04时，参数不敏感[12]。最后采用SCE-UA算法[15]对敏感参数进行优选得到最优参数值，结果见表1。

表1 参数敏感性指数及其率定值Tab.1 Parameter sensitivity index and rating parameters

注：*表示敏感，□表示不敏感。

4 结果分析

由表1中参数的敏感性指数可知，以流量为目标，敏感参数为：土壤与细小裂隙(基质流)面积所占比例(α)、慢速、快速水箱入渗补给常数(ks、kf)、快速水箱出流系数(hF)，以及山坡单元蒸发折算系数(kc)、管道区降水折扣系数(a)。以地下水埋深为目标，敏感参数为：山坡单元的α、a、ks以及洼地单元的hS。采用多目标函数，敏感参数为：山坡和洼地单元α、ks、kf、hS以及山坡单元kc、a。由此可见，无论是以流量或地下水动态变化的单目标，还是反映两者变化的多目标，影响产流的慢速、快速水箱入渗补给参数(ks、kf)以及基质流面积所占比例(α)为敏感性参数；基质区平均张力水容量(wm)、蓄水容量曲线指数(b)均不敏感。

本文对多目标函数敏感的参数进行优选，得到最优参数值(见表1)。山坡、洼地单元之间差异大的参数包括影响水量平衡的蒸发折算系数(kc)、平均张力水容量(wm)以及影响出流过程的慢速、快速水箱出流系数(ηS、ηF)。相比于山坡，洼地土层厚，岩石裸露率较低，基质区面积(α)大，张力水容量(wm)大；由于地形平缓，出流系数(ηF、ηS)小，慢速与快速水箱交换弱(ke小)。山坡管道区最大蓄量VM达70.75 mm，即次降雨量满足VM时产生地表径流，这与该流域坡面集水区实测地表径流形成的累计降雨深基本一致[16]，说明参数率定结果反映坡地、洼地水文地质特征。

流量模拟过程如图3，率定期NSEQ为0.92；实测、模拟径流深分别为300.7 mm、334.0 mm，误差为11.1%；场次洪水过程涨水、退水过程以及峰现时间基本吻合，其中，洪峰最大的洪水过程峰现时间推迟1h，两场大洪水洪峰模拟值偏小，误差分别为-5.6%、-7.7%。验证期NSEQ为0.91；实测、模拟径流深分别为175.6 mm、198.0 mm，误差为12.7%；场次洪水基本吻合，峰现时间基本一致，其中洪峰误差较大的一场洪水误差为10.4%。各场次洪水峰现时间与洪峰误差都在许可误差[6]范围以内。根据模拟结果(图3)，暴雨-洪峰模拟误差相对较大(2017/06/12、2017/07/09)。分析其原因，快速流系统入渗过程随降雨强度变化而变化，尤其在暴雨期间，落水洞集中补给地下管道作用增强，而模型对快速流系统入渗过程进行了概化，未考虑入渗随降雨特征的变化。此外，地下管道汇水过程也随流量变化而变化，模型未考虑汇流过程也是造成误差的重要原因之一。

利用洼地地下水埋深资料对模拟结果进行验证(图4)，率定期NSED为0.81，平均埋深误差为0.01 m；验证期NSED为0.83，平均埋深误差为0.02 m。模拟与实测过程基本吻合。

图3 模拟与实测流量过程与洪水过程Fig.3 Simulated and measured discharge process and flood process

图4 地下水埋深模拟与实测过程线及两者相关关系Fig.4 Simulated and measured groundwater depth process and correlation

山坡和洼地单元蒸散发量以及产流量模拟结果统计值如表2。山坡林地、灌丛等覆被好，实际蒸散发量大，径流深小；但山坡面积大，计算期内产水量占总水量的76.0%。流域地表径流、快速径流、慢速径流占总水量比例分别为3.8%、71.8%、24.4%，地表径流和快速径流量之和远比慢速径流量大，反映了喀斯特地区裂隙管道快速流对陡涨、陡落流量过程的控制作用。

表2 流域水文过程模拟结果Tab.2 Simulation results of Watershed hydrological process

5 结 语

本文针对峰丛喀斯特流域山坡、洼地地貌特征，构建了反映导水介质慢速、快速水流入渗补给和蓄泄特征以及山坡-洼地水力联系的半分布式水文模型，通过典型流域实测数据对模型进行验证。研究结果表明：

模型可以较好地模拟喀斯特地区降雨径流过程、洼地地下水动态变化以及不同地貌单元对流域出口流量的贡献，参数率定值能够反映山坡与洼地地形和水文地质特征对产流过程的影响；山坡单元是峰丛—洼地喀斯特地貌区主要水源，约占总产流的3/4；快速径流约占总产流的2/3，是喀斯特山区陡涨、陡落流量过程主要控制因素。本研究为我国南方喀斯特地区水资源利用和洪涝形成过程预测提供定量分析方法。

模型还有许多不足，模型应用于西南喀斯特小流域，忽略了流域的汇流过程；同时，没有考虑洼地中农业生产活动对产流和地下水过程的影响，在今后的研究中将加以改进。

□