在森林生态系统中, 降水量被分为林冠截留量、穿透降雨量和树干茎流量三部分, 其中林冠截留量在该系统的水文循环和水量平衡中占有重要地位, 约占降水量的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模型的流域模拟中意义重大。

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

1 研究区概况

研究区位于乌鲁木齐县境内的乌鲁木齐河流域(86°45′—87°56′E, 43°00′—44°07′N), 山区流域面积1070 km², 属于中型流域, 见图 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。

2.2 SWAT林冠截留模型

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

(1)

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

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

当R′_day≤can_day-R_INT(i)时：

(2)

当R′_day＞can_day-R_INT(i)时:

(3)

式中, R_INT(i)为初始林冠蓄水量(mm), R_INT(f)为最终林冠蓄水量(mm), R′_day为当日单次降雨事件降雨量(mm), R_day为穿透雨量(mm)。

2.3 Gash模型

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

当林冠截留量未达到饱和时(R′_day＜P′_G):

(4)

当林冠截留量达到饱和时(R′_day＞P′_G):

(5)

式中, p为自由穿透雨系数；p_t为树干茎流系数；R为平均降雨强度(mm/h)；E为饱和林冠的平均蒸发速率(mm/h)；P′_G为使林冠达到饱和的降雨量(mm)。

使林冠达到饱和的降雨量P′_G的计算公式如下：

(6)

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

2.4 Gash模型参数的确定

饱和林冠的平均蒸发速率E(mm/h), 可通过实测数据中降雨量R′_day与林冠截留量I之间的回归关系确定。Gash模型假设在单次降雨事件期间平均蒸发速率E(mm/h)与平均降雨强度R(mm/h)的比值是恒定的^[9], 因此可以通过回归关系的斜率来确定E。但该方法会高估饱和林冠蒸发速率, 所以使用SWAT模型内嵌Penman-Monteith公式输出结果作为饱和林冠的平均蒸发速率, 最终确定E=0.15 mm/h。本研究区平均降雨强度为1.54 mm/h。

林冠枝叶持水能力S(mm), 可以通过穿透雨量和降雨量R′_day的关系方程来确定^[22], x轴上的截距即为林冠枝叶持水能力S, 得出S=2.29 mm(图 2)。

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

树干茎流系数p_t, 可以根据树干茎流量与降雨量R′_day之间的关系方程来确定^[12], 方程的斜率即为p_t, 结合前人研究结果表明天山中段山区树木树皮干裂, 对树干茎流产生很大的吸收滞流作用, 从而较少有树干茎流出现^[24], 所以此处取值0.01。根据上述各参数代入公式(6)计算研究区内主要树种雪岭云杉的P′_G为3.09 mm。

2.5 对SWAT模型子程序的修改

为了对林冠截留过程进行改进, 本研究根据SWAT官网下载的Arcswat源代码Rev.681版本, 使用Inter Visual Fortran对其进行编译和调试后形成可执行文件。SWAT模型由311个子程序构成, 其中每个子程序文件都会对单个或多个过程进行建模。需要修改的部分集中在模型的两个子程序文件中, 在“modparm.f”子程序文件中添加使林冠达到饱和的降雨量P′_G为新的变量, 将公式(4)—(6)在“conopyint.f”文件中进行编译, 带入Gash模型所需要的参数来完成林冠截留物理过程的描述。各子程序文件编译完成后, 主程序“main.f”读取所有输入数据并且调用链接的子程序文件从而进行后续的建模过程。全部子程序运行成功并且生成64位发布模式可执行文件, 将此可执行文件改名“swatUser.exe”放入“TxtInOut“文件中等待调用即可。

2.6 模型校准和验证

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

(7)

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

(8)

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

(9)

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

(10)

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

(11)

(12)

(13)

式中, Q_i^obs为第i个观测径流值(m³/s), Q_i^sim为第i个模拟值(m³/s), Q^omean为观测流量平均值(m³/s), Q^smean为模拟流量平均值(m³/s), σ_sim和σ_obs是模拟值和观测值的标准差(m³/s)。这些指标的可接受范围一般是NSE≥0.6, R²＞0.6, |PBIAS|≤20%。通过泰勒图的形式以RMSD、标准差和皮尔逊相关系数三个指标来量化两种不同SWAT模型观测值与模拟值流量间的一致性。

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 m³/s和6.55—6.86 m³/s之间; NSE值分别在0.63—0.85和0.58—0.82之间；R²值分别在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 m³/s和4.68 m³/s。使用SWAT模型时, RMSD值分别为3.49和7.80 m³/s, 说明与SWAT-Gash模型相比SWAT模型在干旱半干旱区的径流模拟效果欠佳。各点到原点之间的距离表示各组模拟值与观测值的标准差。在校准期SWAT模型和SWAT-Gash模型的标准差分别为7.59和7.33, 验证期标准差分别为6.45和6.54。SWAT-Gash模型在验证期标准差较大, 可能是因为敏感性参数更多, 但总体上该模型仍有更好的表现。相关系数表示出山口实测径流数据和模型模拟值之间的关系, 表明地表径流在时间序列上的相似性。SWAT-Gash模型在校准期与验证期的相关系数分别为0.93和0.81, 高于SWAT模型的0.91和0.77。相较于SWAT模型, 改进模型在校准期和验证期模拟高流量的能力更加突出。

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

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

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模型。