流域水文模型在面源污染模拟与管控中的应用研究进展

谢 晖,邱嘉丽,董建玮,高田田, 4,赖锡军

1 中国科学院南京地理与湖泊研究所,中国科学院流域地理学重点实验室,南京 210008 2 中国科学院地理科学与资源研究所,中国科学院陆地水循环及地表过程重点实验室,北京 100101 3 南京师范大学海洋科学与工程学院,南京 210023 4 南京信息工程大学应用气象学院,南京 210044

面源,或称之为非点源,是以面状分布并排放污染物进而造成水体污染的发生源。随着点源污染控制水平的不断提高,面源污染成为全球水体污染的重要污染来源[1—2]。按照污染类型,面源污染主要分为城市面源和农业面源污染。城市面源污染伴随着城市化的快速发展,不透水下垫面比例增加,暴雨径流及携带的面源污染物超出城市排水系统的消纳和处理能力,对城市水环境带来冲击。近年来雨污分流改造和海绵城市建设不断推进,城市面源污染逐渐得到缓减和控制。我国作为粗放型农业生产大国,对化肥农药的依赖度极高,不合理的水肥管理导致的农业面源污染问题尤为突出,威胁着水环境健康、饮用水和粮食安全[3]。第二次全国污染源普查数据显示,农业源排放的化学需氧量、总氮和总磷分别占水污染物排放总量的49.8%、46.5%和67.2%,表明农业面源污染仍然是造成我国水污染的重要原因。我国高度重视农业面源污染防治工作,自2015年以来,国家相关部委相继出台了《关于打好农业面源污染防治攻坚战的实施意见》、《重点流域农业面源污染综合治理示范工程建设规划(2016—2020年)》等一系列政策性措施,习近平总书记也强调要“以钉钉子精神推进农业面源污染防治”。2021年3月生态环境部和农业农村部联合印发了《农业面源污染治理与监督指导实施方案(试行)》,明确了“十四五”至2035年重点区域农业面源污染防治的工作目标和主要任务,对面源污染的监测和评估提出了更高的要求。

面源污染的发生和迁移包含若干动态过程,如降雨径流、土壤侵蚀、污染物累积与冲刷、土壤溶质淋溶等,具有随机性、复杂性、不确定性及时滞性等特点,获取一定精度和密度的监测数据需要高强度且复杂的野外监测,现有监测网络在点位分布和技术手段等方面还不足以准确评估面源污染对水体污染的贡献,分布式和半分布式流域水文模型逐渐盛行。目前国内外应用于面源污染模拟的此类模型主要有SWAT模型(Soil and Water Assessment Tool)[4]、SWMM模型(Storm Water Management Model)[5]、AnnAGNPS模型(Annualized Agricultural Non-Point Source Pollution Model)[6]、GWLF模型(Generalized Watershed Loading Function)[7]、SPARROW模型(SPAtially Referenced Regressions On Watershed attributes)[8]、MIKE SHE模型[9]和HSPF模型(Hydrological Simulation Program-FORTRAN)[10—11]等。这些模型既有相通之处,又具备各自的独特性,用户可以根据研究区特色、模拟尺度、数据完备性、研究目的等多方面进行模型的适用性选择。本文针对HSPF模型进行综述。HSPF模型是美国环保署重点推荐的用于模拟流域水文及水质的综合模型。它起源于20世纪60年代的斯坦福水文模型(Stanford Watershed Model),随后水质模块被逐渐引入,并最终由Robert Carl Johanson于1980年正式提出[11]。发展至今,HSPF模型的开发和维护主要在美国环保署的资助下进行,现已更新至12.5 版本,并且被集成至综合模拟平台BASINS(Better Assessment Science Integrating Point and Non-point sources)。该平台集成多种模型、数据处理及分析工具,实现了HSPF模型的可视化配置、运行及结果分析。HSPF模型现阶段的主要开发方RESPEC公司强调了该模型的两个特点：“comprehensive(综合)”和“flexible(灵活)”,其综合性体现在模型在连续过程、动态事件以及稳态条件下均可模拟流域水文、河道水力以及水质的复合过程,模拟范围涵盖了面源污染的产生排放、迁移转化及水体输运等主要过程；灵活性体现在模型可根据实际情况调整模拟步长、计算单元的划分方式、陆域与水体的承接关系、污染物迁移的模拟方法等,以应对面源污染在不同尺度需求下的模拟。这两个重要的特点突出了HSPF模型在面源污染模拟方面的优势,使其在全球范围内得以广泛运用。

近年来HSPF模型逐渐被用于我国流域水环境过程模拟,应用过程中也面临着模型机理和参数本地化、模型构建精细化、模型结果不确定性较大、落地应用困难等方面的挑战。本文围绕HSPF模型在面源污染模拟与管控中的研究进展进行了综述。首先简介该模型的结构功能及基本原理,综述其在变化环境下的模拟方法和成果;其次,对应对参数识别、不确定性分析、措施效果评估和总量控制的思路和方法进行了总结;最后,分析了新形势下该模型的延伸拓展,提出了未来HSPF模型可进一步深化的发展动向。

1 HSPF模型结构及功能、基本原理

1.1 结构及功能

HSPF模型的运行主要依赖两个文件,即用户控制输入(User Control Input,UCI)文件和流域数据管理(Watershed Data Management,WDM)文件。前者是ASCII文本文件,以模块化形式记载流域属性和模型配置。WDM文件是一种非格式化、32位架构的二进制文件,它存储HSPF模型所必须的气象驱动数据及模拟结果产生的时间序列数据。在最新版的BASINS 4.5下,这两个文件的生成和编辑可通过SARA Time series Utility工具和winHSPF软件实现,HSPF模型的搭建环境及流程如图1所示。

图1 基于BASINS 4.5的HSPF模型搭建环境及流程Fig.1 Environment and process to develop HSPF model based on BASINS 4.5BASINS: 点源和面源污染评估系统 Better Assessment Science Integrating Point and Non-point Sources；WDM: 流域数据管理 Watershed Data Management；HSPF: 流域水文模型HSPF Hydrological Simulation Program-FORTRAN；UCI: 用户控制输入 User Control Input；DEM: 数字高程模型 Digital Elevation Model

对于一个目标流域,HSPF模型会将其地理空间概化为3类模拟单元：透水陆面(PERLND)、非透水陆面(IMPLND)和水体(RCHRES),每个子流域中则是由不同数量和种类的3类模拟单元组成,子流域间再由拓扑关系相联。HSPF模型对模拟单元划分时可考虑将气象驱动、土地利用和土壤类型等环境因子作为空间异质性的划分准则,并将不同因子映射于每个模拟单元,形成水文响应单元。这种划分方法体现了该模型半分布式的架构,有利于表征面源污染输出的时空分异特征。

HSPF模型功能特点可很好地应对面源污染的模拟。从研究区特征来说,由于HSPF模型对陆面进行了透水和非透水划分,因此可对农业流域、城市流域或混合复杂的用地特征流域进行面源污染模拟。当然,模型对非透水陆面的模拟过程中考虑的物理过程相比于专门针对城市面源污染模拟的SWMM模型来说较为简单,不能体现管网汇流这一城市下垫面特色,但HSPF模型也为模型使用者提供了一种简易的城市面源污染模拟方法。从模拟的时间尺度上看,由于HSPF模型受连续气象数据的驱动,因此可模拟面源污染的长时间连续过程。单次降雨事件带来的面源污染输出同样受到业界的关注,由于HSPF模型模拟时间步长可在小时至分钟调整,输出时间序列可细至每小时,因此也可用于模拟单次降雨事件下的面源污染完整过程。从模拟对象上看,HSPF模型成熟的水文模拟框架能以较高的精度模拟径流和演算河道洪水。悬浮泥沙在陆面的分离搬迁以及在河道中的动力输运过程也在HSPF模型中得以刻画。常见的面源污染物如氮磷的累积冲刷、形态间转化、各相间分配、作物和藻类吸收等基本过程均在模型中有所考虑,因此HSPF模型对面源氮磷污染模拟具有很好的适用性。此外,HSPF模型还有专门模拟农药流失和扩散的PEST模块。对于其他化学物质,HSPF模型依据物质的特性(如保守性、衰减性、吸附性等),可在模块中新增污染物并设定相应的参数进行模拟,已有研究利用HSPF模型模拟过大肠杆菌、汞、有机碳的面源污染过程并取得了较好的效果[12—14]。

1.2 基本原理

径流：径流由PERLND模块和IMPLND模块分别计算,前者是HSPF模型中的核心模块,涉及到的物理过程最多。透水陆面的剖面自上而下分层为地表蓄积、壤中流蓄积、上土壤层蓄积、下土壤层蓄积和地下水蓄积。降水在经冠层截留后落在地表。HSPF模型沿用斯坦福模型的概化方法,将落地雨的潜在去向划分为：形成地表滞蓄增量、形成壤中流滞蓄增量或直接下渗,且入渗量随地块面积呈线性变化。具体的下渗容量分配概化如图2所示。

图2 HSPF模型下渗容量分配(改自Bicknell等[10])Fig.2 Infiltration distribution concept in HSPF model (modified from Bicknell et al. [10])

落地雨经以上分配后,各土壤蓄积进行动态的水量平衡计算。其中从地表蓄积中产生的地表径流以Chezy-Manning方程计算；壤中流按线性水库法计算；补给河道流量的地下水部分计为有效地下水,由地下水蓄积量、地下水蓄泄系数、地下水退水系数、地下水水力坡度构建的非线性方程计算。非透水陆面不考虑土壤分层,产流只包含地表径流,采用的算法与透水陆面中一致。河道和湖库中的径流演算在RCHRES模块中进行,但认为水体具有单向流动性,单层且完全混合。在以上的假定条件下,RCHRES采用水力函数表FTABLE结合插值法模拟水深、容积和流量的关系。

泥沙：HSPF模型考虑了粒径的影响,将泥沙细化为沙粒、粉粒和粘粒三种组分。在透水陆面上对泥沙输移量的模拟过程包括泥沙与土壤基质的吸附、分离,以及分离泥沙和土壤基质的冲刷,冲刷量的计算依赖于地表径流量、冲蚀系数和指数。非透水陆面泥沙输移量的模拟仅考虑堆积和冲刷过程。河道中泥沙首先被区分为粘性泥沙(粉粒和粘粒)和非粘性泥沙(沙粒)。对于粘性泥沙,用户须设置悬浮泥沙的沉降临界剪切力以及床面泥沙的起动临界剪切力。当水流剪切力介于沉降临界剪切力与起动临界剪切力之间时,陆面进入水体的粘性泥沙以及河道中的悬浮粘性泥沙将通过平流过程输运；当水流剪切力大于起动临界剪切力时,河床上的粘性泥沙将被冲起并输运；当水流剪切力小于沉降临界剪切力时,粘性泥沙将沉降至河床。对于非粘性泥沙,模型将首先计算水流对非粘性泥沙的输运能力。当该输运能力超出陆面进入水体的与河道中的悬浮非粘性泥沙之和时,河床上的非粘性泥沙将被冲起并输运；反之,非粘性泥沙将沉降至河床。

污染物：污染物在陆面的迁移过程可通过两种方式进行模拟,一种是基于PQUAL和IQUAL子模块分别在透水单元和非透水单元中进行模拟,模拟方法主要基于经验方程。吸附态污染物的输移通量计算由泥沙冲刷量与冲刷系数的乘积获得,冲刷能力系数是需要率定的关键参数。溶解态污染物的传输路径有地表径流、壤中流和基流,地表径流中的污染物输移通量由地表污染物蓄积量、地表径流量以及地表径流对污染物的冲刷速率共同决定。壤中流和基流携带的污染物质量计算受壤中流和基流中污染物浓度影响,且污染物浓度值可设置为每月变化。另一种是通过NITR、PHOS和PEST子模块,以质量守恒为基础分别考虑氮、磷和农药在土壤蓄积中的转化、淋溶以及迁移等复杂的迁移转化过程。前者允许用户设置多达10种污染物,粪大肠菌群主要通过这种方式模拟。降水发生之前,粪大肠菌群在地表的累积量取决于细菌的积累速率,而最大累积量取决于土地利用类型；随着降水的发生,地表径流的冲刷效应使累积的细菌进入河道,粪大肠菌群负荷与径流量密切相关,而壤中流或地下水的贡献主要根据这些传输介质中菌群的恒定浓度计算。后者是更为复杂的模拟方式,仅能够模拟透水陆面中的氮、磷和农药三种污染物的迁移转化过程。这三种污染物对透水陆面的来源输入主要来自大气干湿沉降和施肥。NITR子模块中考虑的氮的形态包括氨态氮、硝态氮、活性有机氮和非活性有机氮,模拟主要过程包括氮与泥沙间的吸附与解吸、硝化与反硝化、矿化与固定、固氮与挥发、作物吸收与氮回收等。PHOS子模块用于模拟磷酸盐和有机磷在土壤蓄积中储存量的变化,主要模拟过程有磷酸盐与泥沙的吸附与解吸、矿化与固化、作物吸收。PEST子模块用于模拟溶解态、吸附态和晶体态农药在土壤中的吸附、解吸以及降解过程。以上过程采用一阶反应动力学方程、FREUNDLICH方程和温度修正的ARRHENIUS方程等进行表达。

水体的污染物输移过程也可以采用简单或复杂两种模拟方式。简单的模拟方式主要基于通用的GQUAL子模块进行模拟,其中不涉及生物化学过程、溶解态污染物考虑其降解过程、吸附态污染物考虑随泥沙的沉降、悬浮、降解以及吸附和解吸过程等。复杂的模拟方式是为了满足需要考虑氮磷在水体中的生物化学转化作用的用户需求,主要基于RQUAL子模块,该子模块首先模拟溶解氧和生化需氧量的平衡过程,进而模拟无机氮磷的平衡过程,以及氮磷与浮游生物相关的过程,最终计算污染物经随流输移从RCHRES的出口汇出量。

2 HSPF模型在变化环境下的面源污染模拟

2.1 气候变化

气候变化在改变流域水文循环过程的同时,也影响着流域污染物的物理化学过程、迁移转化过程及水体对污染物的稀释能力。气候因子,如降雨、温度是流域面源污染的主要驱动力,因此预测未来气候情景下流域面源污染的发展趋势对未来流域水环境的管理和规划具有重要意义。对气候变化情景模拟的主要思路是将气候预测系统同流域模型相耦合。一般包括以下步骤：流域模型率定与验证；构建气候变化情景；将气候数据输入到模型中,模拟气候情景下的流域水文过程和面源污染输移过程。HSPF模型在评估气候变化对流域水文的影响方面具有很好的适用性。Gizaw等[15]采用HSPF模型评估RCP4.5(Representative Concentration Pathway 4.5)和RCP8.5(Representative Concentration Pathway 8.5)情景下的气候变化对埃塞俄比亚四个主要流域径流量的影响,结果表明气温和降水的升高会导致明显的径流增加。在面源污染方面,张鹏飞[16]通过对密云水库流域内四个气象站1961—2000年40年的气象资料分析,采用统计分析和线性回归的方法预测流域气候变化趋势,并分别从任意情景设置和历史资料两方面确定了 25种气候情景和3种水文情景,之后将不同情景下的气象数据输入至建立的HSPF模型, 模拟密云水库流域不同气候情景下水文与面源污染的变化趋势。Yi等[17]采用HSPF模型评估了RCP4.5和RCP8.5情景下未来2020s、2050s和2080s内单位面积面源总氮和总磷排放强度的变化。Kim等[18]通过全球气候模式(GCMs)预测未来气候情景,结合美国地质调查局的土地利用预测数据,采用HSPF模型分析了气候变化和城市化对博伊西河流域水文和水环境的影响,结果表明气候变化与城市化的综合影响可导致面源污染输出呈现显著的季节变异性。

2007年,美国环保署在BASINS系统中集成了气候评估工具CAT(Climate Assessment Tool),允许用户采用CAT创建的气候变化情景作为HSPF模型的输入,为气候变化对面源污染的影响研究以及水资源管理提供了便利。CAT的主要功能是依据任意基准期的温度和降水变化规律预测未来气候的可能变化情景,基准期气候数据包括历史实测数据与气候模式输出数据[19]。Zhou等[20]采用HSPF-CAT耦合模型评估了未来降水变化对尖峰岭雨林流域水文过程的影响。Taner等[21]在四个GCM的降尺度输出数据的基础上,运用HSPF-CAT耦合模型评估了气候变化对径流量和面源污染的影响,结果表明极端降水事件发生频率的增加提高了洪峰量和营养盐流失量。总之,CAT为气候变化情景的创建与修正提供了便利,通过BASINS平台可较为方便地与HSPF模型耦合,使HSPF-CAT成为评价气候变化对面源污染影响的有效工具。

气候变化对面源污染的影响不仅表现在平均态气候变化的影响上,更高强度的扰动事件对面源污染的产输量有着更为严重的影响,这种极端气候事件包括洪涝、干旱及极端气温等。目前运用HSPF模型评估极端气候事件对面源污染的影响研究相对较少。Ouyang等[22—23]采用HSPF-CAT耦合模型评估了未来可能的气温升高强度(增加1°C和2°C)与极端降水事件(雨强增加10%和20%)对美国亚祖河流域泥沙、硝酸盐和磷酸盐负荷的协同影响。Qiu等[24]结合土壤含水量的变化特征,采用HSPF模型评价了密云水库上游流域干旱与极端降水条件对面源污染的协同作用,结果表明,长期干旱后的极端降水事件可能导致累积的沉积物、营养物、粪大肠菌群与其他污染物大量流失,污染物浓度严重超出水质标准,对下游水体的水环境健康造成威胁。这项研究突出了HSPF模型在评估极端气候对面源污染研究中的适用性,也指出了适应极端气候的流域管理体系的重要性。

2.2 土地利用/覆被变化

土地利用/覆被变化能够带来水文条件、下垫面特征、人类活动强度的空间差异,进而影响污染负荷产生和传输的时空分异。HSPF 模型可被用于模拟历史和未来土地利用/覆被变化情景的影响,也能用于模拟土地利用/覆被与气候变化双重驱动下的流域水文与面源污染输出特征,并表现出较好的模拟精度[25—28]。Bello等[29]采用HSPF模型评估了植被覆盖变化对穆尔河流域水质的影响,结果发现BOD、硝态氮和磷酸盐对植被覆盖变化非常敏感。Bai等[27]通过HSPF模型评估了东江流域土地利用变化下流域面源污染的时空变化特征,发现城市化趋势增加了BOD与总磷的负荷量,城市生活污水与城市雨洪是这两种污染物的主要来源；而流域内总氮主要来源于氮肥的施用,流域耕地面积减少使得氮肥施用量减少,进而导致总氮负荷减少。Liu和Tong[30]采用HSPF模型量化了不同宽度(60 m,90 m和120 m)河岸带情景下径流、硝酸盐及总磷负荷的变化,研究结果指导了河岸带森林和湿地缓冲区的规划方案设计。Mohamad等[31]采用HSPF模型评估了农业用地扩张对面源污染的影响,结果表明农业用地增加60%使得年均总磷负荷和总氮负荷分别增加3.8%和5.3%,农业用地增加2倍使得总氮和总磷负荷量增加约2倍。在太湖西苕溪流域,Luo等[32]首先利用CA-Markov模型生成2020s和2030s的土地利用情景,同时结合GCM的降尺度未来气候情景,协同驱动HSPF模型来模拟流域径流和营养盐输出。

综上,气候变化与土地利用/覆被变化对流域水环境过程有着显著的影响,HSPF模型不仅可有效地模拟单一变化情景下的水文与面源污染过程,还可模拟多种气候变化与未来土地利用变化组合情景对面源污染的影响。

3 模型参数识别与不确定性分析

3.1 参数识别

面源污染过程涉及到水循环、土壤侵蚀、污染物迁移转化等多种物理化学过程,HSPF模型在概化这些过程时所用的参数众多,包括具有明确物理机制的确定型参数(如河道底泥氨氮在厌氧环境下的释放速率等),也包括集总式的概念型参数(如可带走土壤中90%污染物的地表径流速率等),针对前者的取值可通过实测获取,而后者的取值则需要进行参数识别。参数识别是面源污染模型建模过程中的关键步骤,通常包含参数敏感性分析及参数率定。

参数敏感性分析有助于识别对模拟结果具有显著影响的模型参数,在后续参数率定中可重点针对敏感性参数进行调整,可以一定程度上减少模型校验的工作量。此外,由于不同类型参数代表着流域系统内不同子系统过程,被识别的高敏感度参数也有助于甄别关键过程并指导高效削减的管控方案设计。常用于HSPF模型的参数敏感性分析方法包括Morris、LH-OAT、RSA、Sobol′s等[33—36]。在参数调整过程中,需要注意的是模型参数的敏感度可能在不同的水雨情条件和不同时间尺度下发生变化,这就需要考虑参数的动态识别分析[37]。Xie等[38]将HSPF模型、傅立叶振幅敏感性分析、时间变化与多尺度结合的方法结合应用,识别了水沙参数在不同水文情势和不同时间尺度下的敏感度,对关键参数的监测获取及污染削减措施实施提出了指导建议。

针对HSPF模型参数的手动率定,美国地质调查局最早研发了其水文参数的专家决策软件HSPExp[39],之后又研发了该软件的改进版HSPExp+[40],并增加了泥沙和营养盐模块的手动率定验证功能。 Mohamad[41]利用HSPExp+对HSPF模型在马来西亚Johor River上游流域的营养盐模拟进行了率定验证,并以此评估面源污染输出对气候变化和化肥施用情景的响应。现如今越来越多的自动率定技术被应用于HSPF模型的水文模拟研究。早期应用最为广泛的是PEST程序[42],该程序是基于Gauss-Marquardt-Levenberg算法,解决的是模型参数取值寻优过程中的非线性最小二乘问题,由于PEST操作性强且迭代运算次数较少,现在依然在国内流域中使用较多[27,43—45]。遗传算法则是近十年HSPF模型应用研究中被频繁提及的自动率定算法。随着多站点率定、多目标率定的需求,带精英策略的非支配排序遗传算法(Non-dominant sort genetic algorithm-II, NSGA-II)[46]被用于HSPF模型水文参数识别研究中。相较于传统的遗传算法,NSGA-II基于Pareto最优前沿的理论使用非支配占优的排序方法,可求解多目标函数的协同最优解。Xie等[36]针对径流时间序列的特点定义了三种目标函数,利用NSGA-II对连续模拟和降雨场次模拟下的HSPF模型水文参数估值进行了比对。Xie和Lian[47]定义三个水文站流量模拟和观测值的纳什效率系数作为多目标函数,利用NSGA-II对HSPF模型水文参数的最优取值进行了求解。

现阶段关于HSPF模型水质参数自动率定的研究非常有限,事实上随着高频水质自动监测站和嵌套式多站点水质监测的普及,越来越多的河流水质数据可服务于模型的水质参数自动率定。需要注意的是HSPF模型的水质参数中有很多是将多过程集总表达的概念型参数,并不具备共识的参数取值区间,如果局限于采用自动率定算法获得这类水质参数的“最优值”,很有可能出现拟合效果较好但与实际情况偏差较大的情况,未来的研究应在自动率定的过程中充分融合模型使用者的专家经验和先验信息,以保证面源污染模拟结果的可靠性。

3.2 不确定性分析

随着面源污染模型越来越多地服务于水资源管理、污染物总量控制、生态补偿等领域,模型结果的可靠性问题和风险问题备受关注。由于流域内水循环及生物地球化学过程属于降雨-径流、面源污染物迁移转化和河流水质水动力等高度非线性过程的复杂组合,而流域模型对这些机理过程通常概化成一般的数学方程,本质上是对真实过程的一种近似状态的描述,这就会带来模拟结果的不确定性。美国最大日负荷量(Total Maximum Daily Load,TMDL)计划的具体实施中明确规定了安全余量的计算,要求如果基于模拟结果分配负荷则须要考量模拟结果不确定性,以降低流域管理的决策风险[48]。在自然科学领域,不确定性通常被归为随机不确定性和认知不确定性[49],而具体到面源污染模型的不确定性则概括为三类,即参数、输入及结构不确定性[50]。

3.2.1参数不确定性

流域模拟中的参数不确定性最早在流域水文模型中发现和探讨。1992年Beven和Binley[51]认为在模拟水文过程时存在“异参同效”(equifinality)现象,即模型结果受参数组合的综合影响,不同的参数组合可能获得相同或相近的模型响应。为诠释和评估这种参数不确定性的影响,普适似然不确定性评估(Generalized Likelihood Uncertainty Estimation,GLUE)方法被两位学者首次提出,后续被广泛用于流域水文模型中,HSPF模型也得益于该方法在参数不确定性方面得以发展。GLUE方法在应用时需设置模型参数的先验分布和用于评价参数组模拟效果的似然函数,随后对参数组进行抽样并重复模拟,结合似然函数阈值以及参数组的似然权重,最终获得参数组的后验分布以及模拟结果的不确定性区间。HSPF模型水文参数的先验分布通常被假设为均匀分布和三角分布,而纳什效率系数则常被选为筛选模拟结果的似然函数[36,47,52]。水文参数之间复杂的相关关系(例如下层土壤额定含水量与土壤渗透率呈极显著负相关关系,上层土壤额定含水量与地下水退水系数呈极显著正相关关系)也通过GLUE方法的应用被探讨[53]。在水质模拟方面,Roostaee和Deng[54]针对硝酸盐和磷酸根在壤中流、基流和陆面冲刷的相关HSPF模型参数的不确定性进行了探讨。Mishra等[55]则利用GLUE方法针对大肠杆菌的陆面冲刷以及一阶死亡率参数的不确定性进行了评估和量化。Xie等[36]考虑了场次降雨事件下面源污染过程特征,构建了一系列基于污染物在次降雨尺度下传输特征的似然函数并将其融合,用以评估HSPF模型泥沙参数对陆面及河道悬浮泥沙输移模拟的不确定性。这些前人的研究成果为HSPF模型参数不确定性理论和技术提供了宝贵的经验。但作为一种“伪”贝叶斯(pseudoBayesian)方法,学术届对GLUE的一些假设提出了质疑,认为包括非统计意义的广义似然函数的选取和有效性参数组的判断筛选等,会导致参数后验分布的推求过程和结果不具备显著统计特征[56—57]。基于经典贝叶斯理论的马尔科夫链蒙特卡洛(Markov chain Monte Carlo,MCMC)方法[58]是用于面源污染模型参数不确定性分析的另一种方向。MCMC方法的应用同样需要确定模型参数的先验分布,但需定义包含模拟误差概率分布的严格似然函数,构建若干条独立并行的Markov链搜索参数空间,通过参数组的抽样、模拟和对比评价,不断更新样本信息使Markov链收敛,实现具有统计意义的最大后验估计。前有学者对比过GLUE方法和MCMC方法应用于HSPF模型评估大肠杆菌参数的不确定性,结果表明二者推求的参数后验分布与模拟不确定性区间相似,但GLUE方法需要设定较高的阈值,而MCMC方法需要更多的计算资源[59]。类似的对比结果也出现于该方法在其他流域水文模型中的应用研究中[57,60],发现MCMC方法取得的参数后验具有略高的识别度,且可以更准确地将参数不确定性贡献区分于模型总不确定性。后来在高效采样算法的引入下,Markov链的搜索效率得到提升,例如近年来被广泛用于环境模拟的差分进化自适应Metropolis算法(Differential evolution adaptive metropolis,DREAM)以及其改进的DREAMZS算法[60—61],具备多链并行抽样能力,自适应地调整搜索步长和方向[62],使得MCMC方法比GLUE方法可更高效地收敛至后验分布,且收敛结果通常具有更好的参数识别度。但目前集成此类高效算法的MCMC方法尚未见报道应用于基于HSPF模型的面源污染模拟研究中。

总体而言,现阶段HSPF模型的参数不确定性研究多集中于水文模拟,对水质相关参数的研究比较缺乏。由于HSPF模型的水质参数种类众多,涉及陆面、河道、水土中迁移、形态转换等多界面多过程,应借鉴其他面源污染模型应对这种多维度和多模态参数不确定性的先进方法,该方向的研究也将为HSPF模型更好地服务面源污染管控提供理论和技术支撑。

3.2.2输入不确定性

降雨和DEM是面源污染模拟的关键输入数据。针对HSPF模型来说,Wang等[63]认为二者的不确定性会传递至径流、泥沙及氮磷污染模拟并具有一定的扩大效应,该结论也与其他面源污染模型的相关研究结果较为一致[64]。目前降雨数据的获取仍以地面雨量站网观测为主,由于站网布设不合理以及测量误差,会对面源污染模型的模拟结果带来较大的不确定性[65]。增加雨量站点的数量和密度可提高降雨空间异质性的描述,从而提高径流和面源污染模拟的精度,但雨量站点稀疏是很多流域建模时遇到的普遍现象。HSPF模型在BASINS平台中建模时,默认选用的是临近法或泰森多边形法对子流域和雨量站进行匹配,但对于复杂地形流域,可能需要考虑高程、坡度、坡向等要素对降雨插值的影响。Huo等[66]采用信息熵和模糊数学方法,并在考虑了高程数据之后为HSPF模型生成了降雨输入数据,相比常用的泰森多边形法增强了对降雨空间异质性的描述从而提高了径流模拟的可靠性。此外,为分析降雨数据本身的测量随机误差的影响,刘兴坡等[67]以青龙河流域为研究区,区分降雨时间序列的随机成分并建立随机模拟序列驱动了HSPF模型的径流模拟。卫星降水产品由于具有高时空采样频率和覆盖率,也被逐渐用于面源污染模拟的输入。Solakian等[68]基于地面高密度定点雨量站观测数据评价了三种卫星降水产品(TRMM、CRMORPH和PERSIANN)的差异,并通过HSPF模型传递卫星降雨输入不确定性对流量、悬浮泥沙、水温和溶解氧模拟的影响,分析了这种影响在不同模拟空间尺度和季节间的变化。

DEM数据分辨率影响着坡度、坡长、河道坡降等流域关键特征参数的提取,进而改变了相关参数在模型中的输入。DEM数据分辨率带来的不确定性对流域水文及面源污染模拟的重要影响也已在前人的研究中得到共识[69—70]。针对HSPF模型,Roostaee和Deng[71]的研究发现当DEM数据分辨率在3.5 m至100 m区间变化时,地势平缓地区的径流模拟较为敏感,月均流量的相对误差可高达50%；而在坡度陡峭的山区,这种输入不确定性对径流模拟影响非常小。Wang等[63]通过对DEM栅格数据进行重采样获取30 m至480 m的16种分辨率,评估不同面源污染因子对DEM分辨率变化的敏感性。他们的研究结果发现模型对硝酸盐和氨氮的模拟受其不确定性影响最大,且高于对径流模拟的不确定性影响。这表明DEM数据输入的不确定性在径流-面源污染模拟过程中具有传递和扩大效应。

3.2.3结构不确定性

针对模型结构的不确定性分析常用到多模型集合模拟的方法,即加权平均多模型结果以减少单模型由于其对物理过程概化不足造成的不确定性。权重的选取在早期的研究中较为主观且缺乏物理意义,近年来贝叶斯模型平均(Bayesian model averaging, BMA)方法[72]多被应用于面源污染模拟的结构不确定性评价研究,该方法通过似然值代替了权重,具有更高的鲁棒性。Wang等[73]在我国滦河流域评价了HSPF模型模拟氮素污染的结构不确定性,通过集合SWAT模型和GWLF模型并采用BMA方法对多模型结果进行了贝叶斯加权平均。该方法也被用于减少HSPF模型及SWAT模型单独模拟氮磷污染的结构不确定性,通过置信区间的宽度及对观测值的覆盖率等不确定性量化方法提高了模型的解释力[37]。此外,加强模型对机理过程的数学描述方式和方法也是减小模型结构不确定性的一种方向[74],但针对HSPF模型鲜有相关的研究。虽然HSPF模型已经在全球范围被广泛认可,它的模型结构在某些特定现实情景下还有提升的空间。例如在干旱区域非汛期,地表水和地下水的交互作用突出,而HSPF模型尚不能准确表达这种耦合情景。针对水稻田这种特殊的流域地理单元,已有研究对HSPF模型的结构进行了改进,在模型的源代码中引入了圩埂高度和Dirac delta函数以表征水稻田内排水、灌溉和施肥过程,构建了针对水稻灌区的HSPF-paddy模型并模拟了径流、水温、五日生化需氧量、总氮及总磷。该改进模型的模拟结果在田块尺度和流域尺度均得到了验证[75—76]。

4 HSPF模型在面源污染管控中的应用

4.1 BMPs和LIDs效果评估

最佳管理措施(Best management practices, BMPs)和低影响开发(Low impact developments practices,LIDs)是能够有效削减面源污染的管理政策或技术手段[77]。为避免盲目的规划和实施管控方案,需要科学精确地评估BMPs和LIDs的污染物削减效率。面源污染模型是常用的评估工具,研究学者通常采用模型中的措施评估模块或调整措施评估参数的方式,在模型中反映BMPs和LIDs实施的情景并对其效果进行模拟分析。模型模拟相较于实地跟踪监测实验具有较大的灵活性和时效性[78]。BMPs和LIDs的实施单元通常为场地、田间、小区等,通过流域面源污染模型模拟,可将小尺度下的评估结果推广至子流域和流域尺度,也有助于分析管控措施的尺度效应,因此面源污染模型在评估措施削减效果时具有较大的优势[79]。

针对管理型措施,HSPF模型中的SPEC-ACTION模块可用于评估其削减效果。用户在该模块中可设置和调整农业耕作时间、收割时间、化肥施用方案等农业活动的情景,进而模拟免耕、少耕、施肥管理等管理型措施的水质改善效应。针对工程型措施,HSPF模型中内置了BMPRAC模块,用户可自主定义单个工程BMP或LID措施(包括植草沟、过滤带、干湿沉淀池、岸边缓冲带)在流域单元中的数量及污染物去除率[80]。Risal等[81]采用该模块模拟了农田中过滤带及沉淀池在流域尺度下对泥沙、总氮和总磷的削减效率。但BMPRAC模块中直接设置措施效果的处理方式具有主观性,且不能反应措施受降雨、污染源强、地形等多因素影响下削减效果的时空异质性。相比之下,调整措施评估参数则是一种过程型的模拟方法,可以充分结合面源污染模型本身的结构框架,利用模型对污染物迁移转化的描述能力评估措施的污染物削减效果。谢晖[82]分析了四种农业BMPs(免耕、等高种植、梯田和植被过滤带)对面源污染物的削减过程和原理,将措施在流域单元中的数学表达与HSPF模型参数关联,以此评估了三峡库区小流域中四种措施对面源氮磷污染的削减能力。但是这种过程型的措施评估方法与HSPF模型结合的研究较少,相比之下,SWAT模型的在该方向上具有广泛的应用成果。这种差别主要是由于面向用户的易操作性以及应用经验积累的差距。Arabi等[83]首次系统提出了SWAT模型评估农业管控措施的过程型模拟方法；美国德州水资源学院在2011年编撰了基于SWAT模型的管控措施模拟手册[84]；专门用于评估LIDs效果的SWAT-LID模块也被研发[85],这些学者的研究经验和成果都有助于SWAT模型在该领域的推广应用。对于HSPF模型来说,它与SWAT模型有很多共通之处,包括流域地理单元的概化、机理过程的数学表达等等,后续的研究可借鉴SWAT等其他面源污染模型的研究经验,在措施评估模块的开发上有所突破,将有利于HSPF模型高效准确地评估BMPs和LIDs的水质改善效应。

4.2 流域污染物总量控制

流域尺度下的污染物总量控制,通常是指在水质标准的约束下,充分考虑社会经济条件和水体自净能力,控制流域允许排放的最大污染物总量,并对各单元需要削减的污染物负荷进行合理分配[86]。这种水环境管理的理念在美国TMDL计划、欧盟水框架指令中得以发展,并且逐渐在我国多个流域内推广实践。

污染物总量控制须要明确并量化区域污染物排放对水体水质状态的影响,流域面源污染模型则可建立这种驱动-响应关系,并从降雨-径流、面源污染迁移转化、河流水质水动力等多角度分解污染物排放过程。HSPF模型正是当前美国实施TMDL计划的推荐模型,并在Cannon River、Chesapeake Bay、Moosy Creek等多个流域开展了总量控制的可行性研究和实践[13,80,87]。此外,HSPF模型的不确定性对污染总量控制的影响也被研究。Mishra等[55]区分考虑了HSPF模型的随机不确定性和认知不确定性,利用两相蒙特卡洛随机模拟的方法评估了大肠杆菌污染总量控制的情景方案。Ahmadisharaf和Benham[88]结合GLUE方法和HSPF模型,设计了基于水质目标风险的大肠杆菌污染削减情景评价框架,以正向模拟的思路评估了各污染源在不同削减比例下污染物总量超标的概率。另一方面,在一定条件约束下,以反向优化的思路将污染削减量及控制措施分配于不同空间单元成为总量控制领域的研究热点[89—90]。如前所述,HSPF模型具备灵活的面源污染模拟和措施效果评估能力,能够模拟污染削减措施在不同尺度下的环境效益,因此可将其与优化算法或工具进行耦合,实现成本投入-环境效益协同目标下的污染物总量控制。在美国San Antonio River流域内的多个城市型流域,HSPF模型在优化工具SARA Load Reduction Tool和SARA Enhanced BMP Tool的结合下服务了该区域的污染削减和雨洪管理规划[91]。SARA Load Reduction Tool工具可优化分配不同区域的污染负荷削减量,SARA Enhanced BMP Tool可确定控制措施的优化配置(种类、数量和位置等)以确保分配的污染负荷削减量得以实现,HSPF模型则是不同方案制定和反馈迭代过程中的核心过程与纽带。基于反向优化的措施空间配置过程同样受到模型不确定性的影响,Xie等[38]将HSPF模型的参数不确定性引入措施优化配置过程中,在多种水质达标的置信水平下开展了“措施效果评估-不确定性分析-措施优化配置”的一体化研究,研究思路及相关成果服务于减少农村小流域污染物总量控制的决策风险。

5 新形势下HSPF模型的延伸

随着近年来流域水环境综合管理从理论转向实践,国家和地方逐渐推动流域水环境综合管理从信息化向智慧化发展,智慧流域的理念及其信息化平台建设成为重要的方向。地方生态环境部门不仅需要水利、水环境等监测数据的在线掌控,更希望通过仿真模拟实现河湖水质的实时模拟、预报预警和污染源解析等功能,因此对面源污染模型在信息化平台中的集成应用与服务需求愈发强烈[92—93]。HSPF模型在现阶段的操作使用非常依赖BASINS平台,虽然BASINS平台集成了该模型的输入数据处理、模型配置及结果展示等功能,但BASINS的大部分图形用户接口是由Visual Basic 6编译,这与微软Windows现在的主流开发平台Visual.Net并不兼容,且并不能在Windows以外的操作系统(例如类Unix系统)下运行[94]。这些外在因素减弱了HSPF模型自身固有的灵活性以及在信息化平台的可操作性。而对于HSPF模型本身,其最新版本的源代码依然保留了大量Fortran 77语法,内存管理复杂且固化,程序范式严谨但可读性较差,一定程度上限制了该模型的维护、升级和移植。HSPF模型运行的核心输入文件是WDM文件,存储空间有限,且创建和修改需要访问HSPF模型链接库中的许多子程序,用户通常要借助BASINS平台中的组件如WDMUtil、SARA Timeseries Utility等特定工具进行编辑。以上这些HSPF模型固化的特征与现代化环境模拟要求的可复制性和开放性相悖[95],也正是HSPF模型亟需发展的方向。

在流域管理信息化平台建设和现代化环境模拟的背景下, HSPF模型在不同学者和技术人员的共同努力下得到了延伸和发展。PyHSPF[94]是最早的典范,PyHSPF是以Python语言为核心的开源跨平台软件套装,它沿袭了HSPF version 12.2中水和水质模拟的函数库,利用Python联接了函数库、子流域属性和驱动数据,并构建了包含数据预处理、参数自动率定、批量模拟等功能的Python模块和类,为提升HSPF模型模拟效率[96]、模块功能自主研发[97]以及应用编程接口开发[95]提供了多种便利。此外,HSP2(Hydrologic Simulation Program-Python)模型作为HSPF模型的Python版本被研发,该模型使用Python重编译了HSPF模型的核心功能模块,使用跨平台数据存储文件HDF5存储和读写模型的输入和输出数据,并利用交互式Web应用Jupyter Notebook和JupyterLab 维护HSP2的开发、升级和协作[98]。HSP2模型的主要优势在于其模型结构和源代码脱离了特定操作系统和WDM文件的束缚,增强了子流域、土地利用等流域属性和模型参数等的动态表达,从而有利于用户设计多种变化环境的驱动情景,且模型代码可与CPU/GPU异构式并行加速环境兼容。虽然HSP2模型尚处于前期开发阶段,已有研究成功利用该模型在东南亚源头小流域模拟了河流大肠杆菌浓度的高频动态变化[99]。

6 HSPF模型的发展动向及展望

纵观HSPF模型的发展历程,该模型自上世纪80年代被研发后经历了蓬勃的发展和积累。现如今HSPF模型是美国环保署重点推荐的流域水质模型,是各州执行TMDL计划的主要工具之一。虽然该模型在国内的起步较晚,但从最早的径流模拟已扩展至集成多种污染物、多过程的水量水质综合模拟,到如今服务于太湖、巢湖、滇池等流域的污染管控和规划工作,可以说HSPF模型在我国国情和流域背景下的应用有了飞速的发展和实质性突破。尽管HSPF模型具备了多种优点,但现阶段依然存在诸多局限性,作为面源污染模拟的重要工具,该模型的后续完善更应该落脚至解决生态环境领域的实际问题,这就更加需要广大科研工作者的投入和坚持。对于未来HSPF模型的发展动向及展望,笔者认为可以着重关注以下几个方面：

(1)模型本地化改进。HSPF模型虽然已在全球范围内得到成功的应用,但是针对特定的应用场景,如果完全依赖现阶段模型本身的结构功能,还不足以支撑我国河湖管理工作。例如化学需氧量是我国地表水环境质量考核的重要指标,但HSPF模型暂未具有直接模拟化学需氧量的功能模块。此外,闸门和泵站是人类活动干扰流域自然系统的常见水工构筑物,在我国平原圩区十分常见,对此HSPF模型的模拟概化步骤并非友好。因此在未来,应该针对特定的应用场景,改进模型对实际过程的数学表达,或开发全新的功能模块、或耦合其他环境模型,提升HSPF模型在我国流域的适用性。

(2)大尺度精细化模拟。HSPF模型已在大尺度流域(例如166, 000 km2的Chespeake Bay)[64]获得了成功应用,但在大尺度模拟的同时还需考虑计算单元的细化精度。在我国河长制的背景下,计算单元的划分应结合行政区划,服务河长制的网格化精细管理,计算单元由粗糙单元向精细单元转变是大趋势。对于HSPF模型而言,大尺度精细化模拟需要提高计算效率,基于Python的多线程模块和并行计算扩展包如PyCUDA和Numba可以与pyHSPF模型和HSP2模型兼容[96],未来可设计并行架构的云计算环境,将有助于HSPF模型在大尺度复杂流域的应用。

(3)人工智能和模型的互馈集合模拟。当下数字水利、数字流域的建设促进了多源观测资料呈几何级增长,这些基础大数据在人工智能的催化下,可帮助解决HSPF模型所需的高精度输入以及多维度参数自适应寻优问题。此外在少资料地区,流域模型的模拟结果可为人工智能算法提供训练样本[100],从而延伸数据驱动模型的研究范式。因此对于HSPF模型而言,未来应充分结合大数据统计与智能算法,打破学科间隔阂,互馈集合应用以提升面源污染的模拟与预测水平。