SWAT模型在天山林区林冠截留过程中的改进应用

董力轩, 常顺利,*, 张毓涛

1 新疆大学资源与环境科学学院绿洲生态教育部重点实验室,乌鲁木齐 830046 2 新疆林科院森林生态研究所,乌鲁木齐 830063

在森林生态系统中,降水量被分为林冠截留量、穿透降雨量和树干茎流量三部分,其中林冠截留量在该系统的水文循环和水量平衡中占有重要地位,约占降水量的20%—40%[1—3]。因此,在应用SWAT模型[4—5]模拟森林生态水文过程时,林冠截留对地表径流的影响不可忽视。

对林冠截留水文过程的模型化主要可以分为三类,即经验模型、半理论模型和理论模型[6—8]。虽然经验模型可以根据实测数据运用统计学方法直接推求林冠截留量和降雨量之间线性或非线性关系,但不能刻画水文物理过程。因运算繁琐,不易求解,如结合光线与林冠分配降雨规律形成的理论模型在实例中运行困难[6]。以Gash模型为代表的半理论模型相较于前两者工作量小,参数易确定且运算求解方便,显现了一定的优势[9]。Gash模型的适用性较强,能够广泛使用在不同气候类型或林分类型的研究中,如海洋性气候森林[9],热带与亚热带森林[10]、干旱半干旱区灌木林[11]、地中海气候森林等[12]。在国内,Gash模型应用在大兴安岭地区落叶松林[13]、热带季雨林[14]、华北油松人工林[15]、贡嘎山暗针叶林[16]和黄土区人工刺槐林[17]等区域并获得了比较合理的模拟效果。SWAT模型与Gash模型在计算平均蒸发速率时都选择了Penman-Monteith公式,这为两种模型的耦合提供了便利条件[4, 9]。因此,将林冠截留的物理过程引入SWAT模型的流域模拟中意义重大。

本文以天山北坡中段雪岭云杉(Piceaschrenkiana)林为研究对象,结合气象观测、林下穿透雨实测数据和野外调查获得模型所需参数进行建模。考虑到研究区潜在蒸发远远大于实际蒸发[18],林冠截留作用明显,符合Gash模型对每次降雨事件发生时初始条件的要求,进行基于物理过程的林冠截留模拟可以提高模型精度,为水资源管理提供更可靠的依据。

1 研究区概况

研究区位于乌鲁木齐县境内的乌鲁木齐河流域(86°45′—87°56′E,43°00′—44°07′N),山区流域面积1070 km2,属于中型流域,见图1。河源位于天山乌鲁木齐河源1号冰川,自南流向东北,出山口位于英雄桥。根据世界土壤数据库(Harmonized World Soil Database,HWSD),土壤类型有简育灰色土、钙积黑钙土、粘化栗钙土、简育栗钙土、粘化钙积土、永冻薄层土、松软薄层土和冰川。研究区属于温带大陆性气候[18],春秋短暂,夏季昼夜温差大,冬季漫长且多雪。日照时间长且日照辐射大,年平均气温为2—3℃,极端最高气温30℃(7月),极端最低气温-38℃(1月)。降水量少,年降水量400—600 mm,年内降水多集中在6—8月；蒸发量大,年蒸发量980—1150 mm[19]；年均相对湿度64%—66%;干燥度1.2—1.6。该区域植被类型是以雪岭云杉纯林为主的温带针叶林,占流域面积的10%左右,研究区内云杉分布广泛且林冠截留作用明显,流域特征符合研究需求,见图1。

2 研究方法

2.1 模型的输入数据

2.1.1气象数据

选择研究流域范围内的大西沟、英雄桥和小渠子三站1993—2012的日降水、日最高温、日最低温、平均风速等气象数据构建天气发生器,整理后输入气象数据库。

2.1.2数字高程数据

由地理空间数据云http://www.gscloud.cn/下载研究区DEM数字高程数据,选择30 m 分辨率数字高程数据并且转为WGS_1984_UTM_Zone_45N地理投影(图1)。

2.1.3土地利用数据

选择对地物分类更好的哨兵数据,经最大似然分类法得到精度为10 m的土地利用数据(图1)。

2.1.4土壤数据

本研究采用世界土壤数据库,裁剪出研究区,根据SAWT内置土壤类型进行重分类和重采样处理以满足模型的精度要求,使用SPAW软件计算建模所需参数建立土壤数据库。

2.1.5林冠截留观测区布设

依托天山森林生态系统定位站于2010年7—8月设立林冠截留观测区(东经87°27′,北纬43°25′)。该观测区内主要植被为雪岭云杉,无灌木分布。在林间空地布设自记数字雨量仪测定总降水量和降雨强度,在林下设置接雨蓬链接集水桶收集穿透雨量后换算为降雨量深度(mm),具体数据见表1。

表1 林冠层对降雨的分配

图1 乌鲁木齐河流域示意图Fig.1 Sketch map of Urumqi River basin

2.2 SWAT林冠截留模型

当使用SCS曲线法计算地表径流时,不会对林冠截留进行具体物理过程的模拟,估计为产流前地表蓄水量和入渗量的20%。当使用Green & Ampt入渗方程计算地表径流和入渗时,必须额外计算林冠截留量[20]。SWAT允许林冠截留的最大水量每天在变化,该方法根据叶面积指数函数来计算其值。公式如下：

(1)

式中,canday为当日林冠截留的最大水量(mm),canmx为树冠充分发育时的林冠截留量(mm),该参数需要使用者手动输入。LAI为该植被当日的叶面积指数,LAImx为该植被的最大叶面积指数。

在SWAT模型中,得到给定日的降雨数据,模型会通过降雨量的大小来判断是否会有穿透雨形成。公式如下：

(2)

(3)

2.3 Gash模型

Gash模型描述的是一系列彼此分离的降雨事件,每个降雨事件都包含林冠加湿、林冠饱和、以及降雨停止后林冠干燥的过程,且假定每次降雨事件之间有足够的时间让林冠完全恢复到降雨前的干燥程度[9, 21]。模型将整个林冠在降雨过程中各个阶段的截留损失相加得到总的林冠截留量。因为SWAT模型对于林冠截留的计算采用了单次降雨事件单独计算的模式,并且根据降雨量的大小分为林冠截留量未饱和与林冠截留量饱和两种情况,故将Gash模型拆解为分段函数以更好的与SWAT模型拟合。单一降雨事件发生时公式如下：

(4)

(5)

(6)

式中,S为林冠枝叶持水能力(mm)。

2.4 Gash模型参数的确定

图2 降雨量与穿透雨量的回归关系Fig.2 The relationship between throughfall and precipitation

自由穿透雨系数p,为林内测得的平均可见天空率[23],经过野外实地调查郁闭度后确定为0.2。

2.5 对SWAT模型子程序的修改

2.6 模型校准和验证

模型的校准和验证均在英雄桥出水口进行,本文通过决定系数R2、效率系数NSE、相对偏差系数PBIAS和均方根误差RMSD来对模拟值和实测值的拟合度做出评价。决定系数R2的表达式[25]为：

(7)

Nash-Sutcliffe模型效率系数(NSE)的表达式[26]为：

(8)

偏差系数(PBIAS)的表达式为：

(9)

均方根误差(RMSD)的表达式为：

(10)

标准差与皮尔逊相关系数的表达式为：

(11)

(12)

(13)

2.7 不确定性分析

SWAT模拟的流量不确定性主要是由于输入数据的偏差,建模者的专业知识不足或改进模型的结构不合理等原因导致。本文中使用了分位数回归方法,通过计算不同分位数的条件分布来进行非参数不确定性分析[27—29]。确定置信区间后,使用P因子和R因子来量化不确定性。P因子代表了95%的预测不确定性(95PPU),是落在95PPU范围内数据的数量与总数据数量的比值。评价校准质量的另一个参数是R因子,表示95PPU范围的平均宽度与观测数据的标准差的比值。在理想情况下,P因子的取值在0—100%之间,R因子介于0—∞,当R<1.5时认为可接受,在P因子为1且R因子为0是模拟值与观测值完全一致。

3 结果与分析

3.1 参数的敏感性分析

SWAT模型和SWAT-Gash模型均使用SWAT-CUP内SUFI- 2算法进行模拟值的敏感性分析。首先根据类似研究区常用的25个参数[30—31]对SWAT模型进行1000次模拟,发现其中13个参数对SWAT模型模拟结果影响明显(表2),再将其进行500次迭代得到最优校准。同理,对SWAT-Gash模型进行参数的敏感性分析后选择以下16个参数。经分析发现两种模型均对SCS径流曲线系数、土壤饱和水力传导系数、土壤可利用有效水、浅层含水层水位、融雪基温和平均坡长等参数体现出更高的敏感性。

3.2 SWAT模型与SWAT-Gash模型在径流模拟中的性能对比

考虑到气象数据的年际变化,设置1993—1996年为模型的预热期,以减少由于气象数据而导致的潜在不确定性。1997—2004为校准期,2005—2012为验证期。如图3和图4是SWAT模型与SWAT-Gash模型在校准期和验证期的径流模拟过程线,SWAT-Gash模型的模拟值优于SWAT模型的模拟值,体现在SWAT-Gash模型更好的重现了高流量与低流量,SWAT模型会低估高流量。在校准期与验证期,SWAT-Gash模型与SWAT模型在出山口位置的月平均径流量分别为8.45—8.50 m3/s和6.55—6.86 m3/s之间;NSE值分别在0.63—0.85和0.58—0.82之间；R2值分别在0.65—0.86和0.59—0.83之间；PBIAS值分别在9.7%—11.5%和9.2%—17.1%之间(表3)；表明在该研究区SWAT-Gash模型的适用性高于SWAT模型。

校准期和验证期两种模型的均方根误差(RMSD)、标准差(SD)和相关系数等均在泰勒图中体现(图5)。图中各点与x轴上实测值点之间的径向距离表示RMSD的大小,反映了各模型的模拟能力。SWAT-Gash模型得出的RMSD更低,在校准期和验证期的值分别为3.22 m3/s和4.68 m3/s。使用SWAT模型时,RMSD值分别为3.49和7.80 m3/s,说明与SWAT-Gash模型相比SWAT模型在干旱半干旱区的径流模拟效果欠佳。各点到原点之间的距离表示各组模拟值与观测值的标准差。在校准期SWAT模型和SWAT-Gash模型的标准差分别为7.59和7.33,验证期标准差分别为6.45和6.54。SWAT-Gash模型在验证期标准差较大,可能是因为敏

表2 SWAT模型与SWAT-Gash模型敏感性参数

SWAT：水土评估工具 Soil and water assessment tool;CH_K2：主河道有效渗透系数 Effective hydraulic conductivity in main channel;CN2：径流曲线数 Curve number of moisture condition Ⅱ;SOL_K：饱和渗透系数 Hydraulic conductivity of the soil;SOL_AWC：土壤有效含水率 Available water capacity of soil layer;ALPHA_BF：基流α因子 Baseflow alpha factor;REVAPMN：发生再蒸发的浅层含水层水位阈值 Threshold depth of water in the shallow aquifer to “revap” to occur;GWQMN：发生回归流的浅层含水层水位阈值 Threshold depth of water in the shallow aquifer required for return flow to occur;GW_DELAY：地下水时间延迟 Groundwater delay;SMTMP：融雪基温 Snow melt base temperature;SFTMP：降雪基温 Snowfall temperature;SURLAG：地表径流滞后系数 Surface runoff lag time;TLAPS：气温垂直变率 Temperature laps rate;SLSUBBSN：平均坡长 Average slope length; ALPHA_BNK：河岸调蓄基流α因子 Baseflow alpha factor for blank storage;CH_N2：主河道曼宁系数 Manning′s "n" value for the main channel;GW_REVAP：地下水再蒸发系数 Groundwater “revap” coefficient;TIMP：积雪温度滞后因子 Snow pack temperature lag factor;PLAPS：降水垂直变率 Precipitation laps rate;CH_K1：支流有效渗透系数 Effective hydraulic conductivity in tributary channel alluvium;CH_N1：支流曼宁系数 Manning′s "n" value for the tributary channels;DEP_IMP：土壤不透水层深度 Depth to impervious layer in soil profile;EVPOT：壶穴政法系数 pothole evaporation coefficient；EPCO：植物吸收补偿因子 Plant uptake compensation factor

表3 模型径流模拟评价指标

图3 SWAT模型与SWAT-Gash模型校准期径流模拟预测对比Fig.3 Comparison of runoff simulation prediction between SWAT model and SWAT gash model in calibration period

图4 SWAT模型与SWAT-Gash模型验证期径流模拟预测对比Fig.4 Comparison of runoff simulation prediction between SWAT model and SWAT gash model in validation period

感性参数更多,但总体上该模型仍有更好的表现。相关系数表示出山口实测径流数据和模型模拟值之间的关系,表明地表径流在时间序列上的相似性。SWAT-Gash模型在校准期与验证期的相关系数分别为0.93和0.81,高于SWAT模型的0.91和0.77。相较于SWAT模型,改进模型在校准期和验证期模拟高流量的能力更加突出。

图5 SWAT模型和SWAT-Gash模型在校准期和验证期的泰勒图Fig.5 Taylor diagram of calibration period and validation period for SWAT model and SWAT-Gash model

3.3 基于分位数回归的不确定性分析

本文中定义了95%的置信区间β∈[0.025,0.975],使用两种模型的验证期径流量与观测流量进行基于分位数回归的不确定性分析。如图6所示,在验证期SWAT模型的95%置信区间的平均厚度为13.50 m3/s,SWAT-Gash模型为12.86 m3/s,表明SWAT-Gash模型在本研究区的模拟表现更佳。SWAT-Gash模型的P因子为0.96,即96%的观测值可以被该模型95%置信区间所包括,略高于SWAT模型(0.93)。此外,SWAT-Gash模型的R因子(1.19)低于SWAT模型(1.26),表明前者的不确定度范围相对更窄。

图6 SWAT模型和SWAT-Gash模型月径流模拟的不确定性分布Fig.6 Uncertainty distribution of SWAT model and SWAT-Gash model in monthly runoff simulation

4 讨论与结论

4.1 SWAT模型与SWAT-Gash模型在乌鲁木齐河流域的应用

本文将Gash林冠截留模型与SWAT模型耦合后的模拟结果与原始SWAT模型的模拟结果进行对比分析,经过校准和验证后发现SWAT模型与SWAT-Gash模型在该流域的应用均取得了较好的效果,该结论与前人在乌鲁木齐河流域的建模结果一致[30—32],但SWAT-Gash模型表现出更高的模拟精度。当使用SWAT模型在乌鲁木齐河流域模拟时,校准期径流重现表现较好,验证期表现不佳,造成这种现象的原因可能是模型结构复杂,敏感性参数对模拟结果的影响较大[33—35]。该流域内坡面汇流过程、地下水与河道的补给排泄过程和融雪过程[36—37]是影响山区地表径流的主要因素,仅通过调整描述水文过程的参数无法进一步提升模型在该区域的适用性。与SWAT模型相比,SWAT-Gash模型会高估该流域的低流量,Gash模型输入数据的误差也会影响SWAT-Gash模型地表径流的预测。Gash模型计算时要求参数较多,虽然这些参数的确定方法经过国内外前人研究的验证[9,12,17,22—23],可以得到较好的结果,但王艳萍等人[17]研究表明本文方法得出的参数模拟值偏高。Gash模型与SWAT模型耦合时,对林冠截留物理过程的描述更加详细,结合雪岭云杉林下穿透雨实测数据得出满足Gash模型的参数,使SWAT-Gash模型在该研究区的模拟表现更佳。

4.2 SWAT模型与SWAT-Gash模型的不确定性分析

本文采用了基于分位数的不确定性分析。与SWAT模型相比SWAT-Gash模型在基于分位数回归的不确定性分析中的表现更佳,也表明Gash林冠截留模型对SWAT模型在乌鲁木齐河流域的模拟改进效果明显。通常在SWAT模型的建模中,将许多复杂水文过程用简单的经验系数进行参数化,增加了模型的不确定性[27]。而本研究提出的SWAT-Gash模型中,将林冠截留过程和最大冠层截留量参数替换为基于物理过程的方程,从而提高了模型的适用性。SWAT-Gash模型的改进方法可以提高天山北坡中段中小型流域的水文模拟精度,为该区域水资源管理提供更可靠的依据。此外,本研究的降雨量分配数据代表性有所欠缺,还需要在不同地形条件下收集时间尺度更长的降雨分配数据,使用更多样的方法验证Gash模型所需参数,从而在该流域进行更精细化的建模。

4.3 结论

Gash模型参数易于获取,能够通过降雨量、穿透雨量和林冠截留量等数据推求,可以较好的在天山林区进行林冠截留过程的模拟。本文基于SWAT模型,将Gash模型方程代替经验公式提出新的SWAT-Gash模型。此方法可以提高乌鲁木齐河流域水文模拟的精度。同时使用基于分位数的不确定分析发现改进后模型与出山口处实测地表径流数据的相似性更高,适用性高于原始SWAT模型。