栅格岩溶分布式水文模型

2022-01-26周佳奇张企诺陈新宇晁丽君李致家

水资源保护 2022年1期

张珂，周佳奇，张企诺，陈新宇，晁丽君，姚成，李致家

(1.河海大学水文水资源与水利工程科学国家重点实验室，江苏南京 210098；2.河海大学长江保护与绿色发展研究院，江苏南京 210098; 3.河海大学水文水资源学院，江苏南京 210098；4.中国气象局-河海大学水文气象研究联合实验室，江苏南京 210098)

全球岩溶面积约为220万km2，岩溶区居住人口约10亿人。中国西南部是世界上岩溶发育最典型的地区之一，具有集中连片的岩溶地貌，分布面积约54万km2[1-2]。由于石漠化、植被稀少等原因，岩溶流域的土层厚度、土地覆盖类型等下垫面条件空间异质性大，流域内岩溶裂隙发育程度不同，导致了流域产汇流机制复杂，洪水模拟难度大。岩溶流域中的表层岩溶带以及地下深层溶隙、溶洞等都是地表水入渗、地下水赋存和运移的良好介质，在产汇流过程中起着重要作用，探究地表水-岩溶水-土壤水的转化和产汇流规律十分必要。许波刘等[3]将岩溶裂隙概化为5种，推导了5种裂隙的瞬时单位线公式，改进了新安江模型的汇流；Zhou等[4]将集总式新安江模型和两种基于蓄水库的岩溶模型结合起来，建立了岩溶新安江模型(K-XAJ)，在漓江流域的日径流模拟中取得了较好的效果；Jeannin[5]在Holloch洞穴根据观测到的水位、流量、岩溶管道规模以及管道长度建立了水动力模型，该模型可以很好地模拟以管道流为主的岩溶系统；张志才等[6]根据喀斯特流域多孔介质与裂隙水流特征，改进了分布式水文-植被-土壤模型(DHSVM)，并在贵州普定喀斯特水文试验站取得了较好的模拟效果。随着遥感技术在水文领域的不断应用[7-8]，一些学者通过遥感影像分析了岩溶流域石漠化情况[9]。大部分学者将不同发育程度的岩溶裂隙概化为地下溶隙线性水库来构建集总式岩溶水文模型，而构建的分布式概念性岩溶水文模型主要用于日径流和月径流的模拟，构建的岩溶水力学模型需要较为精细的资料，一般用于实验小流域的水文模拟。此外，还有一些研究将机器学习和智能模型用于水文模拟[10-11]。

岩溶流域下垫面空间异质性较大，其洪水过程的精确模拟难度大。集总式的岩溶水文模型难以刻画流域下垫面条件的复杂性[3-4]；水动力岩溶模型在实验小流域模拟效果好，但是需要精细的数据资料，且对于面积更大的中小河流流域，模型的洪水模拟计算量庞大[5-6,12-13]。分布式和半分布式概念性岩溶模型可以在考虑下垫面异质性的同时简化计算，但是目前分布式岩溶模型大部分被用于日径流和月径流的模拟，用于场次洪水模拟的分布式岩溶模型较少[14-17]。针对以上问题，本文在栅格新安江模型(grid-Xin’anjiang model, GXAJ模型)[18]基础上，加入了表层岩溶带自由水蓄水库和深层地下溶隙线性水库，构建了栅格岩溶分布式水文模型(grid-based karst distributed hydrological model, GKDHM模型)，并检验了该模型在典型岩溶流域(普厅河流域)的应用效果。

1 研究方法

1.1 岩溶地貌影响下的降水径流模拟

本文研究的岩溶地貌对流域水文过程的影响主要包括表层岩溶带的调蓄过程以及深层地下溶隙对地下径流的调蓄作用。通过构建表层岩溶带自由水蓄水库和深层地下溶隙线性水库刻画岩溶地区的产汇流过程。

1.1.1表层岩溶带与深层地下溶隙调蓄

图1为降水经过土壤和表层岩溶带调蓄形成水流的示意图。表土层水分下渗进入表层岩溶带，表层岩溶带水分向下渗透，部分进入深层地下水系统，部分由于表层岩溶带下部垂向渗透性减小而形成侧向水流[19]。

图1 土壤与表层岩溶带水流示意图Fig.1 Schematic diagram of water flow insoil and epikarst

基于流域岩溶发育情况和石漠化情况，采用图2所示的表层岩溶带自由水蓄水库模拟表层岩溶带蓄水作用以及出流过程。通过流域岩溶发育情况判断栅格内是否有岩溶发育，无岩溶发育的栅格采用图2(a)的自由水蓄水库来刻画，对于岩溶发育栅格，存在石漠化(图2(b))和未石漠化(图2(c))两种情况，壤中流(QI)为上层土向下层土下渗过程中由于渗透性变弱而产生的侧向径流[20]，石漠化岩溶区土壤层较薄，故认为石漠化岩溶栅格无壤中流，未石漠化岩溶栅格存在壤中流[4]。根据流域石漠化情况，石漠化栅格采用图2(b)的自由水蓄水库来刻画，未石漠化的栅格采用图2(c)的结构。

图1中表土层的调蓄过程通过图2中蓄水容量为SM的自由水蓄水库模拟，图1中表层岩溶带的调蓄过程通过图2中蓄水容量为SMEK的自由水蓄水库模拟，表土层位于表层岩溶带上部，蓄水容量为SM和SMEK的两个自由水蓄水库为串联关系，水流经表层土自由水蓄水库渗入表层岩溶带自由水蓄水库。图1深层饱和水流带中发育的深层地下溶隙，许多研究将其概化为地下溶隙线性水库[3,21]，本文采用同样的方法构建图2(b)和图2(c)中消退系数为CKG的地下溶隙线性水库，模拟深层地下溶隙对地下径流的调蓄作用。

(a) 非岩溶区 (b) 石漠化岩溶区 (c) 未石漠化岩溶区图2 表层岩溶带自由水蓄水库示意图Fig.2 Schematic diagram of freewater reservoirs of epikarst

a.降水经过植被截留，满足张力水蓄水容量后，剩余的雨量进入土壤自由水蓄水库。在岩溶区域内土壤自由水通过渗透过程补给表层岩溶带自由水蓄水库，渗透过程用格林-安普特公式描述[22-23]：

(1)

式中：f(t)为栅格土壤自由水蓄水库下渗能力，mm/h；K为栅格饱和水力传导度，mm/h；Ψ为栅格湿润锋处土壤吸力，mm；Δθ为栅格土壤饱和含水率与初始含水率之差，%；F(t)为栅格累计渗漏量，mm；t为时间，h。

b.表层岩溶带自由水蓄水库对下渗的土壤自由水进行调蓄。定义表层岩溶带侧向径流量为REKI，表层岩溶带自由水蓄水库对地下水的日出流系数为KEKG，REKI和岩溶区域地下径流RG的计算公式如下：

(2)

式中：KEKI为表层岩溶带自由水蓄水库对表层岩溶带侧向出流的日出流系数；SEK为表层岩溶带自由水蓄水量，mm。

c.栅格内表层岩溶带侧向径流采用线性水库进行汇流计算。对于存在岩溶发育的栅格，引入岩溶地下径流消退系数，反映深层地下溶隙与非岩溶区对地下径流调蓄的差异，计算公式为

QEKI,t=CEKIQEKI,t-1+(1-CEKI)REKIU

(3)

QG,t=CKGQG,t-1+(1-CKG)RGU

(4)

式中：QEKI，t、QEKI，t-1分别为t和t-1时刻表层岩溶带侧向出流，m3/s；QG，t、QG，t-1分别为t和t-1时刻地下径流，m3/s；CEKI为表层岩溶带侧向出流消退系数；CKG为岩溶区地下径流消退系数；U为单位换算系数。

1.1.2表层岩溶带蓄水容量空间分布

20世纪70年代Mangin首先用表层岩溶带代表岩溶包气带上部相对含水比较丰富的部分[24]。随着对表层岩溶带的深入研究，表层岩溶带成为岩溶学的重要概念。覃小群等[25]认为表层岩溶带是地表和地下各种岩溶形态构成的不规则带状强岩溶化层。

表层岩溶带在岩溶发育区域普遍存在。对于表层岩溶带的厚度，众多研究表明，在地形平缓的地带，如洼地、岩溶谷地和岸坡等地带，表层岩溶带发育程度较好；在山顶山脊等陡坡地带，表层岩溶带发育规模较小[2,25-26]。

水的侵蚀性是影响表层岩溶带空间差异的重要因素。地形影响区域水循环过程，进而影响岩溶动力过程，对表层岩溶带发育产生重要影响。山坡土壤水有向洼地汇集的趋势，在地势低凹地区土壤含水量相对较高，增加了水流溶蚀时间，促进了表层岩溶带向深部发育[2]。在岩溶峡谷两侧地带和缓坡部位由于其边坡的风化裂隙较其他地貌部位发育强烈，特别是河流两岸岸壁临空面，大气降雨可迅速进入岩体内，大量水流快速进入岩体具有强大的侵蚀力，使得该区域表层岩溶带发育较深。在山顶山脊及陡坡地带，地形陡峭，降雨入渗条件差，雨水多形成坡面流流失，雨水侵蚀的概率较小、时间较短，这些地段表层岩溶带发育较浅[26]。

地形指数反映的是径流在流域中的累计趋势以及重力使径流顺坡移动的趋势[27]，平缓地区的地形指数大，在山顶山脊等陡坡地带地形指数较小[28]。在岩溶发育区域内建立表层岩溶带与地形指数的线性关系，假设岩溶发育区域内地形指数最大的栅格对应的表层岩溶带蓄水容量最大，地形指数最小的栅格对应的表层岩溶带蓄水容量最小，即

(5)

式中：TImax、TImin分别为岩溶发育区域内栅格地形指数中的最大和最小值；SMEKmax、SMEKmin分别为表层岩溶带自由水蓄水容量最大和最小值，mm；ζa、ζb为线性关系的两个系数。

基于数字高程计算得到的每个栅格的地形指数TI，通过率定SMEKmax和SMEKmin并结合式(5)计算得到ζa与ζb后，基于下式可以计算每个栅格的表层岩溶带自由水蓄水容量SMEK：

SMEK=ζaTI+ζb

(6)

图3 GKDHM模型结构Fig.3 GKDHM structure

1.2 GKDHM模型

许多研究表明岩溶流域产流机制为蓄满产流[29-30]，部分适用于岩溶流域降水径流模拟的水文模型是在新安江模型基础上改进而来的[3-4,21]。本文在GXAJ模型的基础上，通过引入表层岩溶带自由水蓄水库构建了GKDHM模型，其结构见图3。石漠化岩溶区土壤层薄，采用二层蒸散发计算方法模拟蒸散发[20,31]，未石漠化岩溶区和非岩溶区采用三层蒸散发计算方法模拟蒸散发。岩溶发育区域土壤下渗补给表层岩溶带自由水的过程采用格林-安普特公式模拟。由于石漠化区域土壤发育差，故认为石漠化岩溶区无壤中流，岩溶区域深层地下溶隙对汇流过程的影响通过引入考虑岩溶区域地下水消退系数的地下溶隙线性水库进行模拟。

图3中地面径流1、壤中流1和地下径流1分别为非岩溶区域的地面径流、壤中流以及地下径流，对应图2(a)中的地表径流QS、壤中流QI和地下径流QG。图3中的未石漠化土壤壤中流2是岩溶区中未石漠化地区的壤中流，对应图2(c)中的壤中流QI。图3中地面径流2、表层岩溶带侧向出流和地下径流2分别是岩溶区的地面径流、表层岩溶带侧向出流和地下径流，对应图2(b)和(c)中的地表径流QS、表层岩溶带侧向出流QEKI和地下径流QG。

GKDHM模型在GXAJ模型基础上新增了以下参数：用于格林-安普特公式的饱和水力传导度fc、土壤孔隙度n、土壤有效孔隙度θe、湿润锋处土壤吸力Ψ(可以根据土壤性质得到)、表层岩溶带自由水蓄水容量SMEK(取值范围10～50 mm)、表层岩溶带自由水蓄水库对表层岩溶带侧向出流的日出流系数KEKI、表层岩溶带自由水蓄水库对地下水的日出流系数KEKG、表层岩溶带侧向出流消退系数CEKI(取值范围0.930～0.998)、岩溶区地下径流消退系数CGK(取值范围0.920～0.998)、表层岩溶带侧向出流马斯京根法汇流参数KEEK、和表层岩溶带侧向出流马斯京根法汇流参数XEEK(取值范围0～0.5)。参数取值参考了新安江模型参数、GXAJ模型参数以及相关研究成果[4,20-21,31-32]。国内外许多研究表明，表层岩溶带的透水性随着表层岩溶带的深度增大而降低[2,33]，对于反映表层岩溶带透水性的参数KEKI和KEKG，假设KEKI与KEKG之和同表层岩溶带厚度呈线性关系，通过率定KEKI与KEKG之和的最大值与最小值，即可确定KEKI与KEKG之和的空间分布，假设KEKI与KEKG的比例为一个常数，通过率定该常数即可确定流域KEKI和KEKG的空间分布。马斯京根法汇流参数参考了姚成[28]对GXAJ模型马斯京根法汇流参数的研究。未介绍的参数可参考GXAJ模型[27，32，34]。

2 研究区概况与数据资料

2.1 研究区概况

以云南省文山壮族苗族自治州富宁县范围内的普厅河为研究区，普厅河属珠江流域西江水系右江二级支流，发源于富宁县花甲乡戈博村附近，于罗村口汇入右江。普厅河流域(23°30′N～23°50′N，105°20′E～105°40′E)是典型的滇东南岩溶高原与广西盆地间的过渡地带。富宁站多年平均径流量1.66亿m3(1984—2014年)，多年平均降水量 1 088 mm(1984—2014年)。普厅河流域为南亚热带季风气候，水汽主要来自印度洋孟加拉湾的西南暖湿气流和太平洋南海或北部湾的东南暖湿气流，具有明显季节性，年平均气温19.3 ℃，最高月平均气温25.2 ℃，最低月平均气温10.8 ℃。流域内建有清华洞水库，该水库为中型水库，于2006年建成验收。富宁水文站以上控制面积为373 km2，流域内设置4处雨量站(图4)。

图4 普厅河流域概况Fig.4 General situation of Puting River Watershed

2.2 数据资料

a.水文气象资料。普厅河流域流量、降水量、蒸发量等水文气象资料由云南省水文局提供，研究所用资料包括日尺度和小时尺度的站点降水、日尺度的站点蒸发、日尺度和小时尺度的富宁水文站流量资料，资料范围为1984—2014年。普厅河流域于2006年修建了清华洞水库，由于没有水库调度资料，并且流域面积较小，流域出口流量受水库修建以及水库蓄泄洪水的影响大，故选取了1984—2003年的数据用于模型模拟。

b.下垫面资料。表1列出了普厅河流域各类土壤所对应的格林-安普特公式下渗参数值。GKDHM模型所需的新增流域下垫面信息为流域岩溶发育情况与石漠化情况。根据地质云(http://geocloud.cgs.gov.cn/)平台的《西南岩溶区石漠化分布图》和《中国南部及东南亚地区岩溶发育程度图》，对地质图进行配准提取和矢量化得到图5。其余所需的下垫面信息与GXAJ模型相同，所有下垫面信息栅格尺度为 1 km×1 km。

表1 普厅河流域格林-安普特公式下渗参数值Table 1 Infiltration parameters for Green-Ampt infiltration equation of Puting River Watershed

(a) 流域岩溶分布

(b) 流域石漠化情况图5 GKDHM模型新增下垫面信息Fig.5 Newly added underlying surface informationrequired by GKDHM

3 结果与讨论

3.1 模拟结果比较

模型模拟结果采用径流深合格率、洪峰合格率、确定性系数3个指标进行评估(表2)，合格标准参考GB/T 22482—2008《水文情报预报规范》。选取1984—2003年的20场洪水，其中前12场洪水用于率定模型参数，后8场洪水用于验证模拟结果。

在径流深模拟方面，GKDHM模型有较好的模拟效果，率定期合格率75.0%，验证期合格率75.0%，相较于GXAJ模型率定期58.3%和验证期50.0%的合格率有显著提升。GKDHM模型率定期和验证期平均相对误差分别为15.5%和17.5%，对比GXAJ模型的18.6%和25.1%，GKDHM模型20场洪水的平均径流深误差降低了4.92%。

表2 普厅河流域模拟结果对比Table 2 Comparison of simulation results of Puting River Watershed

普厅河流域为山区中小流域，较大的洪峰易导致泥石流、滑坡等灾害的发生，能够准确模拟洪峰十分重要。对于洪峰模拟，GKDHM模型率定期合格率50.0%，平均相对误差25.5%，验证期合格率62.5%，平均相对误差20.5%;GXAJ模型率定期合格率50.0%，平均相对误差26.7%，验证期合格率50%，平均相对误差27.5%。GKDHM模型有效地提升了洪峰模拟合格率，20场洪水平均洪峰误差降低了3.59%。

总体来看，GKDHM模型的模拟效果优于GXAJ模型。率定期和验证期的平均确定性系数，GKDHM模型的结果都高于GXAJ模型的结果。20场洪水中GKDHM模型模拟结果确定性系数大于0.5的场次为13场，GXAJ模型为11场。确定性系数的提升表明GKDHM模型能更好地模拟洪水过程。GKDHM模型分布式的结构能更准确地反映岩溶流域空间异质性较大的下垫面情况对产汇流的影响。

表3 普厅河流域模型参数率定结果Table 3 Calibrated parameters of model for Puting River Watershed

3.2 模型参数分析

图6 表层岩溶带自由水蓄水容量空间分布Fig.6 Spatial distribution of free waterstorage capacity of epikarst

(a) KEKG (b) KEKI图7 表层岩溶带自由水蓄水库参数空间分布Fig.7 Spatial distributions of parameters for freewater reservoirs of epikarst

(a) 1985090709号洪水(GKDHM模型) (b) 1985090709号洪水(GXAJ模型)

根据图5(a)，流域内除一般的岩溶发育区外，还存在岩溶洞穴区域。基于式(4)，通过CKG1和CKG2两个岩溶区地下水消退系数分别反映岩溶发育区和岩溶洞穴区域对地下水的调蓄作用。岩溶洞穴退水能力较强，地下水平均汇集时间短，地下水消退系数较小，最终率定的岩溶地下水消退系数CKG2小于CKG1。非岩溶区域的地下水消退系数为CG。

3.3 模拟流量径流成分分析

图8为GKDHM模型和GXAJ模型对两场典型洪水的模拟结果,由图8(b)(d)可以看出，GXAJ模型模拟流量的径流成分以地表径流和壤中流为主，地下径流在整个洪水过程中所占比例很小，且退水缓慢。本流域有大范围岩溶发育并且石漠化范围较大，流域的径流成分中地下径流应该占据一定的比重，且退水较快，因此GXAJ模型径流成分模拟结果并不理想。由于引入格林-安普特下渗过程和表层岩溶带自由水蓄水库调蓄过程，图8(a)(c)显示，GKDHM模型模拟径流成分以地表径流和地下径流为主，壤中流较少，因此GKDHM模型可以更好地刻画岩溶流域土壤发育差、壤中流较少、地下径流比例大和地下径流退水快的特点。