岔巴沟流域汛期径流模拟及地表产流特征分析

2023-01-09景珂星宋进喜庞国伟

水土保持研究 2023年1期

景珂星, 毛欢, 宋进喜, 黄鹏, 吴琼, 庞国伟

(1.西北大学城市与环境学院, 西安 710127; 2.陕西省地表系统与环境承载力重点实验室,西安 710127; 3.陕西省黄河研究院, 西安 710127)

地表径流作为水文循环中不可或缺的基本组成部分，是流域水文过程演变的重要指标[1]。在全球气候变化的大背景下，干旱半干旱地区短历时暴雨强度和极端强降水日数急剧增加[2]，加上人类活动影响下流域土地覆被的空间差异性使得流域径流事件和地表产流能力形成了高度的时空变异性[3]，造成了极端洪涝事件的频繁发生[4]，人民的生命财产安全受到了极大的威胁。为了准确预测径流流量，明确产流规律，流域水文模拟已成为国内外水文科学研究的重要方法。

中国黄土高原地区特殊的土壤条件以及降水变差大且集中的气候特征造成了严重的水土流失[5]。从20世纪50年代末，国家开始在黄土高原采取了一系列的水土保持措施，如修建淤地坝、梯田和植树造林等，使得水土流失得到有效控制，生态环境明显改善[6]。但21世纪以来，随着大量淤地坝淤满，水保措施抵御高强度降水事件的能力下降[7]，当地居民的生命财产安全受到了严重的威胁。因此探讨不同土地覆被条件下的地表产流特征有助于深入了解黄土高原地区降雨—径流的变化趋势和产流规律，可初步为探索流域产水产沙与水土保持服务的权衡关系奠定基础，为推动黄河流域生态保护和高质量发展提供参考。

目前针对黄土高原径流变化及其影响因素的研究，一类是通过径流小区的降雨—径流试验，以控制变量法分析不同地形和降水条件影响下的产流情况[8-9]；另一类是基于流域尺度的水文模型进行径流预测及变化的归因分析[10-11]。然而，由于黄土高原地区流域水文条件时空异质性导致地表径流及产流能力在不同季节和降雨事件内是不断变化的[12]，小范围的降雨—径流试验对物理机制考虑不足，对径流产生和变化的空间规律描述也有一定局限。而SWAT(Soil and Water Assessment Tool)模型不但考虑物理机制且将一些物理机制难以描述的过程用“经验模型”来刻画，这使其广泛应用于不同环境下径流模拟研究。SWAT模型在黄土高原的径流模拟也有广泛的应用，多位学者利用SWAT模型分别在黄土高原汾河[13]、桑干河[14]、渭河[15]、北洛河流域[16]及整个黄土高原[17]的气候变化及土地利用对径流变化的影响进行定量分析；庞佼等[18]利用SWAT模型对黄土高原安家沟流域的月径流进行模拟；Zhao等[19]在泾河流域通过SUFI_2算法对SWAT模型模拟径流的参数进行不确定性分析。这些研究证明SWAT模型在黄土高原区的适用性较好，但多以SCS径流曲线法作为产流方法，与黄土高原超渗产流[20]的产流方式存在较大的差异，且研究尺度多集中于月尺度和年尺度，日尺度径流模拟效果不佳，基于日径流事件的模拟相对缺乏。因此，本文以黄土高原岔巴沟流域为研究区域，采用与黄土高原超渗产流原理相似的Green-Ampt下渗法驱动SWAT模型进行日尺度的汛期(6—10月)径流模拟，并以地表径流系数为指标对降雨强度和前期土壤含水量对不同土地覆被地表产流的影响进行分析，以期推进黄土高原地区日尺度水文模拟的发展，并为揭示黄土高原不同土地覆被产流机制的变化提供参考。

1 研究方法与数据来源

1.1 研究区概况

岔巴沟流域位于陕西省榆林市子洲县北部(109.5°—110.1°E,37.5°—37.8°N)，是大理河的一级支流(图1)，岔巴沟河全长26.6 km，流域面积205 km2，曹坪水文站控制面积187 km2。河流最上游发源于刘新窑村，年平均径流总量约1×107m3，年平均流量约为0.35 m3/s。岔巴沟流域呈对称树叶状，其海拔范围为900～1 284 m，平均海拔大约1 080 m，流域属沟壑密度大的黄土丘陵沟壑区，土层颗粒较小且较为疏松。流域年平均降水量约为474 mm，降水量年内分配极为不均，约70%降雨集中在7—9月，多为历时短、雨强大的暴雨；降水年际变化较大，极易引起山洪泥石流等自然灾害，如岔巴沟2017年“7·26”特大洪水等[21]。目前政府已投入大量资金进行山洪灾害的防控，但相关研究仍落后于治理工程。

图1 岔巴沟流域概况

1.2 数据来源

1.2.1 水文气象数据本研究所用的径流数据选取位于流域出口的曹坪水文站的汛期日径流数据；降水数据主要选取流域内13个雨量站的汛期降水摘录数据。降雨径流数据来源于《黄河流域水文资料—黄河中游区上段(河口镇—龙门)上下册》，时间范围为2006—2013年。由于流域内没有气象站，根据需要选取流域周边绥德、横山、榆林3个国家气象站的逐日气温、蒸发、风速和相对湿度数据，由中国气象科学数据共享服务网(http:∥data.cma.cn)提供。水文站和雨量站点流域内分布情况见图1。

1.2.2 遥感数据构建SWAT模型需要的遥感数据主要包括流域数字高程模型(DEM)、土地覆被、土壤类型等。本研究采用的流域数字高程模型是30 m分辨率的GDEM数据，由中国科学院地理空间数据云平台(http:∥www.gscloud.cn)提供；土地覆被数据为2010年30 m分辨率的中国土地利用遥感监测数据，由中国科学院资源环境科学数据中心(http:∥www.resdc.cn)提供；土壤数据为1 km格网的HWSD(Harmonized World Soil Database)世界土壤数据库[22]。

1.3 SWAT分布式水文模型

1.3.1 SWAT模型原理及结构 SWAT模型对流域的水文建模主要分为水文循环的陆地阶段和汇流阶段。SWAT可以直接从文件中读取输入数据，如降水、最高和最低气温等，从而在运行时生成模拟器。SWAT模型对水文循环的模拟基于水量平衡方程，其表达式为：

(1)

式中:SWt表示土壤最终含水量(mm);SW0表示第i天土壤初始含水量(mm);t表示时间;Rday表示第i天的降水量(mm);Qsurf表示第i天地表径流(mm);Ea表示第i天的蒸散发量(mm);Wseep表示第i天从土壤剖面进入包气带的水量(mm);Qgw表示第i天基流回归的水量(mm)。

本文以Green-Ampt下渗法代替SCS径流曲线法作为SWAT模型的产流方法，Green-Ampt方程假设地表永远存在多余水分进而预测下渗，模型假设在湿峰以上的土壤完全饱和，在湿峰处含水量急剧变化。Mein等[23]开发了一种方法来确定Green-Ampt方程的下渗蓄水时间，该方法被引入SWAT中作为模拟地表径流的替代方法，需要用户输入时间尺度小于日的降水数据计算径流。Green-Ampt Mein-Larson方程的下渗速率定义为：

(2)

(3)

式中:finf,t为时刻t的下渗速率(mm/h);Ke为有效水力传导率(mm/h);Ψwf为湿润锋土壤基质势(mm),是土壤孔隙率、砂粒百分比和黏土百分比的函数;Δθv为湿润锋土壤体积含水量的变化(mm/mm);Finf,t为时刻t累计下渗量(mm);Ksat为饱和导水率(mm/h)；CN为径流曲线数。

1.3.2 SWAT模型构建构建SWAT模型的前期准备为土壤和气象属性数据库的建立。土壤数据库包括岔巴沟流域各类土壤的物理和化学性质，土壤物理性质中土壤的颗粒组成在HWSD数据库中直接获取，并作为SPAW软件的输入数据计算土壤湿密度(SOL-BD)、土壤层有效持水量(SOL-AWC)和饱和导水率(SOL-K)等；土壤化学性质主要包括pH、电导率(SOL-EC)、有机碳(SOL-CBN)等，由HWSD土壤数据库直接获取。气象数据库包括气象站气象数据以及岔巴沟流域13个雨量站降雨数据，以.txt格式输入模型。

通过基础数据的准备和数据库的构建，借助ArcSWAT 2012版本软件建立岔巴沟流域SWAT模型。首先基于DEM数据生成河网，通过推荐阈值将岔巴沟流域划分为31个子流域，并设置流域总出水口；其次通过叠加土地利用数据、土壤数据和坡度数据，将三者的最小阈值均设置为0%，保留流域所有具有相同土地利用类型、土壤类型和坡度的640个水文响应单元(Hydrologic Response Units)；最后将气象数据导入模型，改变流域产流方法为Green-Ampt下渗法后运行SWAT模型，实现岔巴沟流域2006—2013年水文过程的模拟，模拟结果“output.hru”文件中的时间步长内HRU总降水量(PRECIP，mm)、产生的地表径流量(SURQ_GEN，mm)和初始时段土壤剖面的含水量(SW_INIT，mm)等数据用于后续地表产流特征分析。

1.3.3 SWAT模型敏感性分析和参数率定不同的参数对SWAT模型的模拟结果影响很大，因此在模型率定前需进行参数的敏感性分析，选取合适的参数以提升模拟效率。本研究利用SWAT-CUP中的SUFI-2(sequential uncertainty fitting Version2)作为优化算法进行参数敏感性分析和率定。参数敏感性分析结果采用t统计量(t-Stat)和显著性指标p值(p-value)来衡量，t统计量绝对值越大、显著性指标p值越接近0，参数敏感性越强[24]。表1列出了岔巴沟流域率定使用的23个参数及其详细信息。

表1 参数含义及初始率定范围

模型效率通过确定系数(R2)、Nash-Sutcliffe效率(ENS)和百分比偏差(PBIAS)第3个指标来进行评估。指标的计算公式如下[25]：

(4)

(5)

(6)

在3类指标中，ENS的范围为-∞到1，较高的ENS值是首选；R2的取值范围从0到1，数值越大表示模拟效果越好；PBIAS比较模拟输出与观测数据的平均趋势。一般情况下，当R2，ENS>0.5且PBIAS≤±25%，则认为模型效果较好。

2 结果与分析

2.1 岔巴沟流域SWAT模型敏感性分析结果

为了得到模型相关各参数的敏感性信息，本研究基于SWAT-CUP软件采用SUFI-2算法对模型23个有关参数进行了多次迭代，得到了各参数的敏感性分析结果和参数率定范围，其中敏感性较高的参数有CN2，CANMX，SLSOIL，LAT_TTIME，SOL_AWC和SOL_K等，表2列出了敏感性参数的验证范围、t统计量和显著性指标p值的详细结果。敏感性参数同Gao等[26]在黄土高原泾河流域的敏感性分析结果相近。其中CN2(径流曲线数)是对径流影响最显著的参数，CN2值越大表明流域下垫面不透水性越强，产生的径流量越大，CN2的t统计值为-7.239，p值为0，是最敏感的参数；树冠最大储量CANMX影响降雨的树冠截留量，值越大径流量越小；SOSOIL为侧向地下水流坡长，数值大小决定流域坡面特征对HRU径流计算的影响程度；LAT_TTIME为横向流动时间，受土壤水力性质影响很大，影响流域汇流速度；SOL_AWC和SOL_K是与土壤属性相关的参数，其中SOL_AWC为土壤层有效含水量，值越大土壤可容纳水分越多，径流量越小；SOL_K为土壤饱和水力传导率，值越大土壤渗透性越强，径流量越小。

表2 参数敏感性分析结果

2.2 岔巴沟流域SWAT模型模拟结果

在岔巴沟流域SWAT模型径流模拟过程中，将2006—2007年为模型预热期，以2008—2010年汛期为模型率定期，2011—2013年汛期为模型验证期。以R2，ENS，PBIAS第3个指标对模型的适用性进行评价，表3是流域总出水口曹坪水文站汛期日尺度径流率定期和验证期的模拟评价结果，结果表明曹坪水文站的汛期流量模拟值与实测值吻合较好，率定期和验证期的ENS系数分别为0.76，0.74，R2分别为0.78，0.75，模拟结果较好。PBIAS的计算结果表明模型在率定期将流量低估了11.3%，而在验证期将流量高估了17.5%，这可能是由于验证期特别是2013年的实测日径流量变化很大，急起急落现象突出造成的。图2为曹坪水文站率定期、验证期的模拟结果，由图可以看出SWAT模型可以较好的模拟流域汛期径流的变化趋势。

图2 曹坪站汛期径流实测值与模拟值比较

模型参数的不确定性采用P-factor和R-factor这2个指标来衡量，其中P-factor表示95 PPU区间包含实测数据的百分数，R-factor定义了95 PPU的频带平均宽度除以相应测量变量的标准偏差。理论上，P-factor的范围是0～100%，R-factor的范围是0～∞，P-factor接近1，R-factor接近0是完全接近于实测数据的模拟[27]。表3中模型率定期和验证期P-factor分别为0.52，0.51，表明50%以上的实测数据在95 PPU置信区间范围之内，而R-factor分别为0.32，0.25，表明模型率定的准确性可以接受，模型参数的不确定性较小。总体而言，以Green-Ampt下渗法驱动的SWAT模型在岔巴沟流域表现出较好的模拟精度，模型率定期和验证期ENS系数均大于0.70，R2均达大于0.75，PBIAS介于±25%，参数不确定性较小，可以进一步分析流域不同土地覆被下的地表产流特征。

表3 模型模拟结果评价

2.3 岔巴沟流域地表产流特征分析

为进一步理解流域地表径流的形成机制和变化规律，引入可以综合降雨特征、前期土壤含水量和地表特征的变量—地表径流系数作为评价指标[28]，在汛期径流模拟的基础上，基于率定好的模型，以次降雨事件内水文响应单元平均降雨量大于25 mm且平均地表产流在2.5 mm以上为条件筛选出四次降雨—径流事件，结合640个水文响应单元累计得到2 560个不同降雨量—产流量的分析对象，用于分析不同土地覆被下的地表产流特征。

2.3.1 降雨强度对不同土地覆被产流特征的影响降雨强度是黄土高原丘陵沟壑区产流产沙的主要影响因子[8]。根据最大1 h降雨强度等级和土地利用类型将2 560个分析对象的地表产流量和径流系数进行分类统计得到不同最大雨强等级[29]下耕地、草地、林地的地表产流量和径流系数(图3)。从整体趋势来看，随着最大雨强的增加，3种土地覆被的地表产流量和径流系数都逐渐上升，尤其当最大雨强从8.1～15.9 mm/h提高至16 mm/h以上时，3种土地覆被下地表产流量和径流系数均显著增大；不论在何种最大雨强区间下，3种土地覆被地表产流量以及径流系数的关系均为耕地>草地>林地，且随着最大雨强的增加三者间的差距也逐渐加大。

从不同土地覆被来看，地表产流随不同土地利用类型的变化差异较明显，最大雨强在2.5 mm/h以下时耕地、草地和林地均没有产生地表径流，耕地和草地在最大雨强达到2.6 mm/h后开始产生地表径流，最大雨强在8 mm/h以上林地才开始出现产流；随着最大雨强的增加耕地平均地表产流量从0 mm逐渐提升至0.91 mm，4.14 mm和12.83 mm，径流系数在最大雨强在16 mm/h以上时达到0.28，是受雨强影响最显著的土地利用类型，其次为草地，林地变化幅度最小，这是由于森林可以通过减小雨滴动能、拦截雨量以及改变土壤结构和地表结皮来削弱雨强对地表产流的影响[30]。综合来看，流域地表径流系数随最大雨强的增加呈显著上升趋势，且最大雨强达到16 mm/h以上时地表径流系数显著增加。在同等级雨强下，耕地将降雨转化为径流的比例最大，其次为草地和林地。

图3 不同降雨强度下土地覆被地表产流能力对比

2.3.2 前期土壤含水量对不同土地覆被产流特征的影响前期土壤含水是指暴雨事件发生前流域表层土壤的含水量，在温带气候下的许多研究表明前期的土壤水分条件是产生径流的一个重要因素[31]。为了进一步分析不同土地利用类型下前期土壤含水对产流的影响，选取产流能力最强的降雨强度区间，四次事件累计出现498个最大雨强在16 mm/h以上的分析对象，其中耕地214个、草地197个、林地87个，得到其前期土壤含水量与径流系数的拟合关系见图4。从图4A可以看出498个分析对象的初始土壤含水与地表径流系数拟合R2为0.35，说明前期土壤含水量的大小可以揭示大雨强下部分地表径流系数的变化，地表径流系数随前期土壤含水量的增加呈上升趋势。

将分析对象按照土地覆被类型细分后拟合关系进一步提高(图4B，C，D)，其中林地前期土壤含水量与地表径流系数的拟合关系最好，R2为0.47，其次为耕地和草地；从拟合方程的斜率来看，3种土地覆被下地表径流系数对前期土壤含水量的响应存在差异，其中耕地和草地的斜率分别为0.049，0.050，大于林地的0.044，说明耕地和草地地表径流系数对前期土壤含水量的变化敏感程度更高，即随着前期土壤含水量的增加地表径流系数的变化幅度比林地更大。这可能是由于森林涵蓄土壤水分以及截留降水的作用减少了地表径流的产生，造成其对前期土壤含水量变化敏感程度较低。

图4 不同土地覆被前期土壤含水量与地表径流系数拟合关系

3 讨论

以Green-Ampt下渗法为产流方法的SWAT模型在参数率定后达到了较好的模拟效果，但部分学者的研究结果认为降雨数据的时间分辨率对Green-Ampt下渗法的径流模拟效果影响很大[32]，由于研究区内的雨量站均为汛期站，仅获得了5—9月2 h步长的降雨摘录数据，对降雨的强度输入有所削弱，加上流域周边气象站点分布较少使得气象因子的输入精度有一定偏差，对模型产流有一定的影响。DEM数据的空间分辨率对坡度、河网及子流域划分等建模过程有显著影响，进而影响到流域径流模拟精度，有研究表明20～150 m分辨率可以基本满足流域模型模拟精度的要求[33]，但岔巴沟流域属于典型的黄土丘陵沟壑区，地形破碎，沟壑纵横，在陕北黄土高原的研究发现利用DEM提取的平均坡度随数据分辨率的降低呈线性下降的态势[34]，因此需要更高精度的DEM数据来识别研究区的坡度和沟壑密度。本研究由于未能获取更精细的DEM数据，故使用30 m分辨率的DEM数据进行水文分析，对径流模拟的精度可能会造成一定影响。在之后的研究中将利用无人机实地获取分辨率更高的DEM数据来降低地形对径流模拟的影响。