GWLF模型的改进及其在新安江流域的应用研究

2020-12-21齐作达亢戈霖王玉秋

水资源与水工程学报 2020年4期

齐作达，亢戈霖，王玉秋

(1.交通运输部天津水运工程科学研究所，港口水工建筑技术国家工程实验室，工程泥沙交通行业重点实验室，天津 300456；2.南开大学环境科学与工程学院天津市跨介质复合污染环境治理技术重点实验室，天津 300350)

1 研究背景

当今社会，水污染问题越来越受到人们的重视。近些年来，世界范围内许多国家和政府逐步增加了对于水资源保护的投资。考虑到地表水监测人力、财力和物力的花费，非点源流域模型作为一种关键而有效的决策工具应运而生[1-3]。

非点源污染已经成为当前流域污染的重要组成部分，对于这类污染来源的控制与削减迫在眉睫[4-5]。与点源排放不同的是，非点源污染具有间接性的特征，并且与水文过程密切相关，因此识别起来更加困难和复杂[6-7]。为了能够更加准确而简洁地描述流域中水文循环和污染物归趋，恰当的模型框架和适用性高的方程显得尤其关键[8]。到目前为止，人们已经提出了许多能够概化流域水文和污染物循环的方程，如SCS-CN径流曲线方程、Green-Ampt方程、Muskingum洪水演进方程等[9-11]。

通用流域污染负荷方程GWLF(generalized watershed loading function)是一个集总参数流域负荷模型，被美国TMDL计划所推荐和采纳，并且在全世界各地得到了广泛应用。模型在考虑一定的水文循环和污染物归趋的基础上，对于数据的需求量较小[12]，同时可以在月尺度上提供较为可靠的模拟结果，能够为管理决策提供极大的便利和支持。近些年来，GWLF不仅在美国、欧洲等地区取得了较好的模拟结果[13-16]，在国内新安江流域、天津市于桥水库上游流域、滦河流域等地区也有了诸多成功的应用案例[17-20]。

本次研究在水文框架改进的基础上，进一步为GWLF模型添加了河道传输过程中泥沙的沉积、悬浮过程以及污染物的降解过程。在新安江流域的应用结果表明，改进版GWLF模型与原始GWLF模型相比，能够更加准确地模拟流域内流量、泥沙量、溶解性氮负荷量和溶解性磷负荷量[21]。

2 GWLF原始模型及其改进方法

2.1 原始GWLF模型

GWLF模型是一个集总参数流域模型，模型需要输入逐日的气温和降水数据，然后在日尺度上对流域的水文平衡进行连续模拟，最终累计到月尺度并输出水量结果，泥沙和营养盐则基于月尺度模拟并逐月输出模拟结果。模型不考虑河道传输过程，并假设当月流域内产生的污染物全部传输到流域出口。水文产流计算采用了SCS-CN径流曲线方程，泥沙计算公式采用了USLE通用土壤流失方程，营养盐计算使用了经验系数方程。更详细的模型理论描述可以参考GWLF理论参考手册。

2.2 GWLF模型的改进

2.2.1 GWLF模型的改进方法 GWLF模型的改进包括以下几个方面：(1)将整个流域划分为若干子流域，每个子流域内地表和地下的水文和营养盐独立计算，并汇入子流域对应的河道内；(2)模型的水文、泥沙、营养盐基于日尺度计算；(3)添加河道传输过程，河道来水包括上游汇入、本子流域地表和地下水汇入以及点源的排放，河道的水文循环使用曼宁方程和变量存储法进行描述，采用巴格诺德简化方程表征泥沙在河道中的悬浮和沉积过程，营养盐衰减使用与滞留时间相关的一阶指数方程进行描述。改进后的模型框架示意图见图1。

图1 改进后的模型框架示意图

2.2.2 GWLF改进模型的水文、泥沙、营养盐相关计算

(1)河道日平均流量计算

河道的日平均流量计算公式如下[22]：

(1)

将河道假设为河岸坡度为2∶1的梯形(见图1)，河道中的水深可以用如下方程计算：

(2)

式中：H为河道水深，m，其下标j为迭代次数，H0= 0；Bi为河段i的河底宽度，m；ni为河段i的曼宁系数;Slopei为河段i的河岸坡度。

河段内的平均流速和传输时间可以使用如下方程表示：

(3)

TravelTimei,t=lengthi·1000/(Vi,t·3600)

(4)

式中：Vi,t为第t天河段i的平均流速，m/s；TravelTimei,t为第t天河段i的传输时间，h；lengthi为河段i的长度，km。

存储系数SCi,t使用如下方程计算：

SCi,t=min(1,48/(2·TravelTimei,t+24))

(5)

最终河段的流出水量为FlowOuti,t(m3)：

(6)

(2)河道泥沙浓度计算

河道中的泥沙浓度计算公式如下：

(7)

式中：SedCi,t为第t天河段i的泥沙浓度，t/m3；SedSti,t为第t天河段i进水后的泥沙总量，t。

河道中的峰值流速计算公式如下：

velocityPki,t=pvaf·Vi,t

(8)

式中：velocityPki,t为第t天河段i的峰值流速，m/s；pvaf为峰值径流调整系数。

河段最大可传输泥沙浓度计算公式如下：

SedCmaxi,t=lineari·velocityPki,texponenti

(9)

式中：SedCmaxi,t为第t天河段i的最大可传输泥沙浓度，t/m3；lineari和exponenti分别为线性系数和指数系数。

河段的泥沙沉积或悬浮量计算公式如下：

(10)

式中：SedDepi,t为第t天河段i的泥沙沉积量或泥沙悬浮量，t。

(3)河道中的营养盐计算。

河道中的营养盐衰减过程描述方程如下：

(11)

3 研究区域概况及数据来源

3.1 研究区域概况

本次研究选取新安江屯溪监测站上游流域作为研究区域，研究区域概况如图2所示。所选研究区面积约为2 683 km2，是中国第一个跨省生态补偿示范区，包括率水和横江两条主要支流(图1(a))，长度分别为97 km和55 km，主河道总长385 km，下游汇入千岛湖，为杭州市主要饮用水源地。流域被分为13个子流域，共有14个雨量监测站(图1(b))。流域海拔在96～1 622 m之间(图1(c))，气候以亚热带湿润季风气候为主，1997-2013年期间年平均气温为17.2℃，年平均降水量为1 741 mm，大部分降水集中在3-8月。研究区域内的主要土地利用类型及面积占比为林地81%、耕地14.9%、草地3%(图1(d))。

图2 研究区域概况图

3.2 数据来源

本次研究的14个降水站逐日数据和屯溪监测站的逐日流量、逐月泥沙数据来自《中华人民共和国水文年鉴》，逐日的气温数据来自中国气象数据网，数字高程数据(DEM：30 m×30 m)来自地理空间数据云，土地利用数据来自国家第二次土地利用调查数据，逐月的溶解性氮和溶解性磷监测数据由黄山市环境监测站提供。模拟时间范围为1997-2013年，其中1997-1999年的数据用于模型预热。具体数据的格式和详细内容参见表1。

表1 数据类型及来源

4 结果与分析

4.1 模型校准方法与结果讨论

模拟值与实测值之间的比对评价，使用了3个国际公认的统计学指标，分别为纳什系数(NSE)、决定系数(R2)和观测标准偏差率(RSR)，其计算公式如下：

(12)

(13)

(14)

模型的校准方法采用了基于贝叶斯理论的广义似然不确定性估计法(GLUE)，按照流量、泥沙、溶解性氮、溶解性磷的顺序依次校准相关参数，综合3项统计学指标选出最优参数组，校准后得到的参数数值如表2所示。

由表2可以看出，原始GWLF模型与改进模型校准后的各参数数值存在差异，但是差别较小，主要参数如径流曲线系数、通用土壤流失方程系数、地表径流氮、磷浓度和地下水氮、磷浓度差值保持在5%范围以内，两个模型的参数结果均在合理范围之内，校准后的参数组可以用于模拟屯溪流域的流量、泥沙和营养盐负荷量。

表2 原始GWLF与改进后的模型参数校准结果

原始GWLF模型与改进后模型的各项模拟结果的统计学指标如表3所示。

表3 模拟结果的各项统计学指标

由表3可知，改进模型的流量模拟统计学指标优于原始GWLF模型，NES、R2值均大于原始GWLF模型，尽管模拟逐日尺度流量结果不是GWLF模型原有的功能，但是改进模型的NES依然在0.9左右，表明改进模型在屯溪流域能够提供较为精确的逐日流量模拟结果；原始模型泥沙的NES和R2均小于0.8，而改进后模型的NES和R2均大于0.9，同时改进模型的RSR值小于原始GWLF模型，表明改进模型在计算泥沙方面有所优化；溶解性氮的模拟结果统计学指标差异与泥沙相似，改进模型要明显优于原始GWLF模型，但是溶解性磷的模拟结果统计学指标方面，改进模型与原始GWLF模型的差别较小，除RSR以外，另外两项指标改进模型仅高出不到5%。

4.2 流量模拟结果

2000-2013年逐月流量实测值与模型模拟结果对比如图3所示。从图3整体上看，两个模型的模拟计算结果令人满意，几乎再现了实测流量的变化趋势，且模拟值与实测值的峰值保持一致。此外，改进模型在枯水期的模拟效果要明显好于原始GWLF模型，原始GWLF模型倾向于低估枯水期的流量，这种差异在月尺度的模拟结果中并不是十分显著。但是在日尺度上，二者的流量模拟结果差异较大，改进后的模型极大地提高了模拟的准确性。

图3 2000-2013年逐月流量实测值与模型模拟结果对比

由于2000-2013年14 a的日尺度模拟结果时间序列图清晰度较差，所以选择以2005年为例，展示两个模型之间逐日流量模拟结果的详细差异，如图4所示。由图4可见，无降雨时期，两个模型均无法完全重现流量的变化趋势，但是原始GWLF模型过分地低估了流量，几乎接近于零，并且过于平缓，对于较少降水的响应较差；而在大多数情况下，改进后的模型模拟结果与实测数据非常接近，并且基本捕捉了流量的变化趋势。这种差异的重要原因是改进模型增加了河道过程，并引入了子流域概念。原始GWLF模型还存在高估暴雨事件流量峰值的现象，同时也会高估随后的流量，这主要是因为没有河道过程，缺少河道存水，当日的子流域产生的径流全部转化为流域出口的排水量。

图4 2005年逐日流量实测值与模型模拟结果对比

4.3 泥沙量模拟结果

2000-2013年逐月泥沙量实测值与模拟结果对比如图5所示。由图5可看出，两个模型对泥沙量的模拟效果较好，基本捕捉到了泥沙产量的波动性。但是原始GWLF模型对峰值模拟准确性较低，往往极大地高估或低估了极值，改进后的模型模拟值更加接近于实测值，其对峰值模拟结果有了明显的改善。结合表3中的统计学指标可知，改进后GWLF模型对于泥沙的模拟效果要明显优于原始GWLF模型。这种差距说明模型水文循环结构改进优化了泥沙模拟结果，增加河道泥沙悬浮沉积过程，使得改进后模型能够对河道泥沙的携带能力做出相应调整，模拟效果更好，适应性更强。

图5 2000-2013年逐月泥沙量实测值与模拟结果对比

4.4 溶解性氮模拟结果

2000-2013年逐月溶解性氮负荷量实测值与模拟结果对比如图6所示。

图6表明，与实测数据相比，原始GWLF模型和改进模型的模拟结果重现性较好，均能很好地再现实测数据变化趋势。但原始GWLF模型对极值的模拟效果不甚理想，倾向于高估峰值而低估极小值，在枯水期模型响应较差，对低值的波动性捕捉较差。改进后的模型不仅优化了水文模拟过程，还新增了污染物河道衰减指数，对峰值和低值的模拟更加准确，波动极端性得到改善。同时改进后的模型模拟结果统计指标NES和R2均大于0.9，RSR小于0.3，明显优于原始GWLF模型模拟结果。修正后的模型在一定程度上解决了原始GWLF模型对极值模拟的缺陷问题，模型模拟效果得到进一步优化。

4.5 溶解性磷模拟结果

2000-2013年逐月溶解性磷负荷量实测值与模拟结果对比如图7所示。由图7可见，虽然对溶解性磷的模拟结果与溶解性氮相比不甚理想，但两个模型基本上可以反映溶解性磷的历史变化趋势。原始GWLF模型对峰值模拟往往偏低，平水期模拟偏高，枯水期时间序列过于平缓，无法准确捕捉到观测数据的每一次波动。改进后模型对极值的模拟在一定程度上弥补了原模型的不足，模拟效果得到了一定的改善。这种差异可能来自于模型机理的改进：在水文过程优化的基础上，将原始GWLF模型中城市营养盐负荷为固态的假设改为城市总磷营养盐负荷的1/4为溶解性磷[23]，并添加了污染物河道衰减过程。