基于GWLF模型的于桥水库流域氮磷负荷估算及来源变化解析

李泽利，罗娜，2，苏德岳，高锴*，赵兴华，2，梅鹏蔚，张震

（1.天津市生态环境监测中心，天津 300191；2.天津天滨瑞成环境技术工程有限公司，天津 300191；3.天津市引滦工程于桥水库管理处，天津 301900）

于桥水库位于天津北部蓟州区城东4 km 处，是引滦入津输水工程的重要调蓄水库，同时也是天津市工农业和生活用水的重要水源地[1]，自1983 年引滦入津以来一直为天津市人民生产和生活提供优质水源。然而自2014 年开始，于桥水库富营养程度加剧[2]，2015—2017 年夏秋季库区暴发大规模蓝藻水华，严重威胁了天津市城市供水安全。氮、磷作为浮游植物赖以生长的重要营养物质，参与光能转化代谢过程，是导致湖库富营养化最为重要的两个因子[3-4]。由于氮、磷的输入与流域内人类活动密切相关，若要控制湖库富营养化和蓝藻水华，首先要明确进入到水体中氮、磷营养物质的来源及负荷[5]。于桥水库的汇水主要来自引滦调水，沙河、黎河、淋河3 条入库河流和库周30 余条沟渠及大气降水等。研究表明，引滦源头潘家口、大黑汀水库水质恶化，以及流域氮磷面源污染加重，是于桥水库水质恶化、水生态平衡遭到破坏的重要原因[6-8]。为了保护和改善于桥水库水质，“水十条”实施以来，津冀两省（市）政府采取多项举措来加强流域的保护力度[9]：①2016 年双方签订《关于引滦入津上下游横向生态补偿的协议》[7]，天津市根据跨界断面水质达标情况拨付补偿金；②2017 年河北省对潘家口、大黑汀水库养鱼网箱进行全面拆除[10]；③天津市在入库沟渠截污、农村污水纳管和农业面源污染防治方面投入了大量资金，流域氮磷营养盐输入情况发生了较大变化。解析于桥水库流域氮磷营养盐来源、输入负荷及变化趋势，对区域水环境管理具有重要意义。

目前国内外在面源污染负荷估算上，机理模型的使用具有主导地位[11]。赵磊[12]总结了国内外应用广泛的GWLF、SWAT、AGNPS 等模型的主要优缺点和适用范围，表明GWLF 模型具有结构简单、参数易于获取等特点，同时该模型的水文模块原理与SWAT、AGNPS 相同，在我国北方滦河流域[13-14]、南方新安江流域[15]都取得了很好的应用效果。本研究首先利用GWLF 模型的水文学模块模拟于桥水库各个子流域的降雨产流量，然后乘以子流域出口断面的氮磷浓度，进而估算氮磷面源污染负荷，最后针对2016—2017 年于桥水库氮磷污染的来源、负荷和变化趋势进行分析，以期为流域水环境管理提供技术支持。

1 材料与方法

1.1 研究区概况

于桥水库正常蓄水面积为86.8 km2，控制流域面积2 060 km2，其中上游河北省境内面积占79%，下游天津市境内面积占21%。流域属于温带半湿润大陆性季风气候，四季分明，多年平均气温12.5 ℃，平均降水量670 mm，主要集中在6—9 月。流域内主要入库河流为沙河、黎河（引滦调水通道）和淋河，此外库周还有30余条入库沟渠，多为季节性河流。

研究区内设有4 个水文站，分别是大黑汀分水闸、前毛庄、水平口和淋河桥站，其中大黑汀分水闸站观测引滦调水流量，前毛庄站观测黎河自产流量与引滦调水流量，淋河桥站观测淋河子流域自产流量，但受上游龙门口水库蓄水影响长期断流，水平口站观测沙河流域自产流量，常年有水流动。研究区内设有4个水质监测断面，分别是隧洞出口、黎河桥、沙河桥和淋河桥，其中隧洞出口为大黑汀水库经引滦暗渠进入于桥水库流域的第一个监测断面，黎河桥、沙河桥和淋河桥分别为黎河、沙河、淋河的跨省界监测断面（图1）。此外，库周还设置38 个沟渠水质观测点[16]，逐月开展水质监测工作。

1.2 研究方法

于桥水库流域氮磷污染来源从空间上可分为引滦调水、沙河、黎河、淋河、库周地表水、库周地下水和库区大气降水7 部分，从管理角度可汇总为引滦调水源、入境河流源（沙河、黎河和淋河）和境内源（库周地表水、地下水和大气降水）3 大部分。用流量乘以氮磷浓度是计算负荷的最直接方法，但由于流域内各条河流的水文站和水质断面的位置在空间上并不全部重合，且沙河、黎河下游和库周均无流量监测数据，无法直接计算。为此，本研究选取气象、水文、下垫面数据比较齐全且受人类干扰较少的沙河水平口子流域，使用GWLF模型模拟产流过程并率定水文参数，然后将水文参数推及到整个流域[17]，模拟每个水质断面子流域产流量，再乘以相应断面的氮磷浓度，进而估算整个流域的氮磷污染负荷。

图1 于桥水库流域河网及子流域示意图Figure 1 Sketch map of river network and subbasins in Yuqiao reservoir watershed

于桥水库流域水量和氮磷污染负荷估算所需的数据见表1。

1.2.1 GWLF模型原理

GWLF 模型是1987 年由美国宾夕法尼亚大学开发的一个半分布式、半经验式的流域氮磷污染负荷模型[18]，该模型对地表径流量的计算基于日气象数据，使用径流曲线数模型（SCS-CN）计算目标流域内不同土地利用类型的径流深，作为表征地表径流量的模拟结果。其中最重要的参数是表征地表产流能力的径流曲线数（Curve number，CN）值，该值与土地利用类型、土壤类型密切相关。模型对地下潜流的计算基于日水量平衡，以浅层饱和区非重力自由水的迁移转化为核心，将浅层饱和区看作一个简单的线型水库进行模拟，通过两个关键参数退水系数r（指浅层饱和区水体横向转移成为地下潜流水的比例）和渗漏系数s（指浅层饱和区水体纵向向下转移进入深层饱和区的比例）调整水量平衡。将逐日的地表径流量与地下水潜流量相加，可得到逐日的河川径流量，即到达流域出口的产流量[15]。

1.2.2 GWLF模型构建

GWLF 模型产流模拟的输入数据主要包括气温、降水、土地利用类型及面积等。首先基于ArcGIS 平台加载研究区的数字高程图（DEM），根据实际进行河网矫正后，以各个水质监测断面为出口分别生成沙河桥、黎河桥、淋河桥子流域和库周子流域（图1），然后与土地利用类型图叠加，分别提取出每个子流域中各个土地利用类型的面积，作为模型的输入数据；气温和降水数据以日为步长作为模型输入数据，水平口水文站实测流量数据以月为步长作为模型校准数据。与产流过程有关的参数，如各类土地利用类型的径流CN值、地下水退水系数r、渗漏系数s和不同月份的蒸发覆盖因子等，先采用模型推荐值输入，然后进行率定。

1.2.3 GWLF模型参数率定

利用沙河水平口水文站月观测流量进行水文参数率定，模拟时间跨度为2006年1月—2017年12月，共144 个月。使用美国康奈尔大学基于GWLF 模型原理开发的ReNuMa 软件[15-16，19]进行参数校准，该软件在Excel文件下运行，利用了Excel强大的规划求解宏，基于月实测数据校准一个或多个参数，使得模拟结果与实测结果尽可能接近（误差平方和最小）[20]。将2006 年的数据作为模型“预热期”，使模型初始条件稳定平衡；2006—2013 年作为模型校准期，校准的参数主要有径流曲线数CN 值、地下水退水系数r和渗漏系数s等；2014—2017 年作为验证期，验证模型参数的可信程度。采用Moriasi 等[21]推荐的统计指标和评价标准判断模型性能和适应性，包括判定系数（R2）、纳氏系数（NSE）、百分比偏差（PBIAS），其评价标准见表2，后两者的计算公式如下：

式中：Qoi为观测值；Qpi为模拟值为多年平均观测值；n为样本个数。

如表2 所示，水平口子流域径流模拟除了验证期的NSE＜0.75 外，在全模拟期的R2、NSE 和PBIAS 分别为0.77、0.76和-3.4%，模型性能总体表现非常好。从图2 的2006—2017 年的降水过程与径流拟合效果图来看，模拟流量和实测流量过程线吻合得较好，GWLF 模型适用于于桥水库流域的径流模拟。使用率定好的模型分别在沙河桥、黎河桥和库周沟渠子流域建模，通过变换土地利用类型和面积等输入数据模拟各个子流域出口的产流量，其中库周沟渠子流域模拟出的地表径流和地下潜流分开参与下一步的负荷量估算。

1.2.4 氮磷负荷估算

表2 基于月尺度径流模拟的模型性能等级评价标准及评价结果Table 2 Evaluation standard of model performance levels and results based on monthly runoff simulation

图2 沙河水平口子流域降水过程与径流拟合效果图Figure 2 Precipitation process and runoff fitting diagram of Shuipingkou subbasin

引滦调水、沙河、黎河、淋河的氮磷负荷分别根据大黑汀分水闸水量、沙河桥模拟水量、黎河桥模拟水量、淋河桥站实测水量乘以隧洞出口、沙河桥、黎河桥、淋河桥断面对应月份的总氮（TN）、总磷（TP）浓度计算得出。2017 年库周地表水氮磷负荷估算已另文发表[16]，2016年采用相同方法进行估算。库周地下水中TN 浓度取自蓟州区2019 年24 个千人以上地下水水源地水质监测结果，平均值为14.3 mg·L-1；TP 浓度取自蓟州区2016—2017 年2 个国家考核地下水井监测结果，平均值为0.04 mg·L-1。根据地下潜流水量模拟结果乘以上述氮磷浓度值，估算出库周地下水输入到于桥水库的氮磷负荷。库区降水输入氮磷负荷采用降水量×水面面积×降水中氮磷浓度来估算，降水中TN 浓度根据2016 年天津市降水离子成分监测结果折算得到：雨水中铵离子（以N 计）浓度年均值为2.62 mg·L-1，硝酸根离子（以N 计）浓度年均值为1.30 mg·L-1，因此TN（以N 计）年均浓度为3.92 mg·L-1；降水中TP 浓度无实测数据，根据太湖流域研究经验[22]，按照0.05 mg·L-1计算。

2 结果与讨论

2.1 氮磷负荷及来源构成

表3 为2016—2017 年于桥水库流域产流量、氮磷负荷估算结果。从表3 可以看出，2016 年于桥水库流域产流量约4.6 亿t，TN 输入负荷2 185.9 t，TP 输入负荷150.2 t；2017年产流量约4.1亿t，TN 输入负荷2 963.1 t，TP 输入负荷39.0 t。徐宁等[23]利用SWAT 模型模拟2010—2013 年于桥水库入库污染负荷，结果显示年均TN 入库负荷为2 900.6 t，TP 负荷为163.2 t，本研究估算的TN、TP 负荷与其在量级上一致，表明GWLF模型适用于于桥水库流域污染负荷研究。

表3 2016—2017年于桥水库流域产流量、氮磷负荷估算结果Table 3 Estimation results of water yield and nitrogen and phosphorus load of Yuqiao reservoir from 2016 to 2017

将2016—2017 年于桥水库流域产流量、氮磷负荷估算结果求平均值后进行来源解析（图3）。从图3可以看出，产流量由大到小依次表现为引滦调水＞沙河＞黎河＞大气降水＞库周地表水＞淋河＞库周地下水，分别占37.3%、18.3%、12.2%、12.0%、9.0%、7.7%和3.5%；TN负荷贡献由大到小依次表现为引滦调水＞沙河＞黎河＞淋河＞库周地表水＞库周地下水＞大气降水，分别占27.9%、25.3%、10.8%、10.7%、8.9%、8.4%和8.0%；TP负荷贡献由大到小依次表现为引滦调水＞库周地表水＞沙河＞黎河＞大气降水＞淋河＞库周地下水，分 别 占66.3%、13.0%、9.2%、7.0%、2.8%、1.1% 和0.6%。

图3 不同来源产流量及氮磷负荷贡献比例Figure 3 Contribution ratio of water yield and nitrogen and phosphorus load from different sources

从管理角度将氮磷负荷估算结果归并为引滦调水源、入境河流源和境内源3 部分，并与天津市水务局2010 年发布的《于桥水库周边水污染源近期治理工程实施方案》有关研究结果[24]进行比较。从表4 可以得出，与2010年相比，2017年入境河流源TN、TP贡献比例分别下降10.2 个百分点和9.2 个百分点，但总体贡献均在30%以上；引滦调水源TN、TP 贡献比例分别上升2.9 个百分点和10.6 个百分点，总体贡献均在30%以上；境内源TN 贡献比例上升7.3 个百分点，TP 贡献比例下降1.5 个百分点。值得注意的是，2016年引滦调水TP 贡献比例高达75.4%，还需进一步开展研究分析。

表4 基于管理角度的于桥水库氮磷负荷来源比较Table 4 Comparison of nitrogen and phosphorus load based on management perspective of Yuqiao reservoir

2.2 氮磷来源变化分析

氮磷负荷受水质和水量共同影响。由于在估算库周地表水、地下水和大气降水等境内源氮磷负荷时，2016年和2017年氮磷浓度数据相同，其负荷量变化主要受年降水量影响，2016 年流域平均降水量为675.6 mm，2017 年为637.7 mm，降幅为5.6%，2017 年境内源氮磷污染负荷均低于2016年。对2016—2017年引滦调水、沙河、黎河、淋河输入的氮磷负荷变化进行分析如下。

2.2.1 引滦调水污染负荷变化

引滦调水水质和水量共同决定调水污染负荷。从图4 可以看出，引滦调水TP 浓度由2016 年初最高的0.74 mg·L-1明显下降到2017 年的0.09 mg·L-1，降幅87.8%。同时，调水水量由2016 年（1 月、5 月、6 月之和）的1.91 亿m3下降到2017 年（6 月、8 月、9 月之和）的1.34 亿m3，下降了29.8%。水量和浓度的共同下降使得TP 负荷由2016 年的113.3 t 下降到2017 年的12.1 t，降幅高达89.3%，这也是表4 中引滦调水贡献比例由75.4%下降到31.1%的主要原因。资料显示2016 年11 月—2017 年5 月，国家、河北省、天津市财政共同出资清理潘家口和大黑汀水库的养鱼网箱，共清理拆解网箱79 687个，拆除后大黑汀水库总磷浓度明显下降[10]。

图4 2016—2017年引滦调水量和污染物浓度变化Figure 4 Comparison of water diversion and pollutant concentration in Luanhe River from 2016 to 2017

然而，引滦调水TN 负荷不降反升。监测结果显示引滦调水TN 浓度由2016 年初的2.18 mg·L-1上升到2017 年的7.39 mg·L-1，增幅239%，这可能是由于调水期8—9 月正值汛期，此时段滦河流域面源污染较重[25-26]。尽管2017 年调水水量较2016 年下降了29.8%，但调水中TN 浓度的升幅较调水水量的降幅变化更大，导致2017年引滦调水TN 输入负荷较2016年上升了123.3%。

2.2.2 入境河流污染负荷变化

沙河、黎河、淋河子流域的面源污染物通过降水产流过程进入于桥水库。由于三个流域具有相似的面源污染特点，现以沙河子流域为例分析TN、TP 负荷变化。沙河子流域2017年产流量较2016年上升了11.4%。从图5 可以看出，2017 年沙河各月TP 浓度均低于2016 年，平均浓度由0.238 mg·L-1下降到0.069 mg·L-1，下降了71.0%，由于TP 浓度的降幅较产流量的升幅变化更大，相乘后使得2017 年TP 负荷较2016年下降了51.9%。2017 年沙河各月TN 浓度与2016年相比互有高低，平均浓度由9.95 mg·L-1下降到9.29 mg·L-1，下降了6.6%，但由于2017 年产流量的升幅较TN 浓度的降幅变化更大，相乘后使得2017 年TN负荷较2016年上升了20.5%。

图5 2016—2017年沙河子流域产流量及氮磷浓度变化Figure 5 Comparison of water yield and TN，TP concentration in Shahe River from 2016 to 2017

2016 年津冀两省（市）政府签订《关于引滦入津上下游横向生态补偿的协议》以来，河北省加大了沙河、黎河、淋河流域污染治理力度，河道水质明显改善，TP 等五项考核指标年均值均满足《地表水环境质量标准》（GB 3838—2002）Ⅲ类水质要求[7]。但由于TN 浓度尚未纳入考核体系，关注程度和治理力度不大，流域内农业种植[27]、畜禽和水产养殖、农村生活污水排放的氮素很容易通过水体流动输入到于桥水库，为藻类大量生长提供充足的营养物质。

3 结论

（1）GWLF 模型对于桥水库流域水文模拟效果较好。2016—2017 年于桥水库流域年均产流量约4.4亿m3，其中引滦调水贡献占比37.3%，3条入境河流贡献38.2%，库周产流和大气降水共占24.5%；年均总氮负荷2 574.6 t，其中引滦调水贡献占比27.9%，3 条入境河流贡献46.8%，库周产流和大气降水共占25.3%；年均总磷负荷94.6 t，其中引滦调水贡献占比66.3%，3 条入境河流贡献17.3%，库周产流和大气降水共占16.4%。

（2）与2016年相比，2017年于桥水库流域总磷负荷减少了74.0%，主要是由于引滦调水和入境河流总磷浓度均明显下降（降幅分别为87.8%和71.0%）。然而，全流域总氮负荷却增加了35.6%，这是由引滦调水总氮浓度大幅上升（增幅239%）以及入境河流产流量增加（增幅11.4%）导致。

（3）于桥水库流域总磷污染防控取得了一定效果，但总氮污染有加重的趋势，需要加强流域非点源氮污染防治力度，同时选择大黑汀水库水质较好的月份调水，以减少氮污染负荷。