1 引言(Introduction)

水质模型已经成为当前水环境规划和管理的重要工具之一, 常见的可以用于模拟水质的模型包括：与农业过程相关的APSIM、DRAINMOD、GLEAMS等；主要应用于流域尺度的SWAT、HSPF、INCA等；针对土地侵蚀过程的USLE、WEPP等；着重关注河道的QUAL系列、WASP、EFDC等(Fu et al., 2019).然而, 能够独立完成河网水力和水质模拟的只有MIKE 11等少数几个模型.MIKE 11作为一款商业软件, 功能齐全、方便准确, 但应用成本较高.

HEC-RAS模型(Hydrologic Engineering Center River Analysis System)可以免费下载使用, 适用于一维恒定流/非恒定流模拟、泥沙迁移、水质分析等, 在水力方面应用广泛(Khattak et al., 2016; Adams et al., 2018), 并且已被证明可以解决河网水质模拟问题(王宇骏等, 2015; 邹丽芬, 2017).有研究利用HEC-RAS模型对河流中的总磷、总氮和化学需氧量等进行了模拟(Haddout et al., 2016; 李国华, 2018; Beck et al., 2019), 也有研究将HEC-RAS应用于最佳管理措施(BMP)的确定、水质相应系数的计算中(Littlejohn et al., 2014; 董飞等, 2016).但目前利用HEC-RAS模型进行水质模拟的研究相对较少, 且多应用于单一河段, 缺少在平原河网地区的应用.

模型参数的取值直接影响水质模拟结果(Zwart et al., 2018).由于水质模拟涉及参数较多, 往往不可能对所有参数进行实测, 因此模型在应用前需对参数进行率定.参数率定分为手动率定和自动估值, 试错法是最常用的手动率定法, 但该方法主观性强、效率较低且具有不确定性.自动估值法可有效避免这些问题, 提高模型求参的效率和准确性, 目前已有研究中提出了采用遗传算法、粒子群算法、近似贝叶斯估计等方法开展模型的参数率定(任婷玉等, 2019).在自动率定时通常借助Matlab软件、Python语言对模型参数进行率定(Leon et al., 2016; Zwart et al., 2018).本研究拟采用Matlab软件对HEC-RAS模型中的水质参数进行率定, 探求一种便捷的平原河网水环境模型建模方法.

2 材料和方法(Materials and methods) 2.1 研究区域概况

研究区域位于浙江省湖州市长兴县中心城区, 主要包括雉城街道、太湖街道和龙山街道、画溪街道、小浦镇的部分区域.该区域人口分布密集, 地势低洼, 河道纵横交错, 池塘星罗棋布, 闸坝、桥梁、排涝站等水工建筑众多(图 1).区域内县级河道包括合溪新港、长兴港、张王塘港、北塘港等8条, 其余均为乡镇级河道, 多数直接或间接汇入合溪新港和长兴港.合溪新港和长兴港主要流向为由西向东, 但当太湖水位过高时, 也会出现倒流或滞流.其他南北向河流的流向基本是由南向北或由北向南分别汇入合溪新港与长兴港.

2.2 研究方法 2.2.1 模型简介

HEC-RAS模型由美国陆军工程兵团水文工程中心研发, 于1995年公开发布1.0版本, 在4.0版本中增加了泥沙运移和水质模拟功能, 目前的5.0版本主要包括恒定流/非恒定流模拟、泥沙输移和水质分析3个模块, 能够应用于水面线计算、河道构筑物影响分析、河道输沙量计算、水质模拟等领域(Hydrologic Engineering Center, 2016), 本研究使用的是HEC-RAS的5.0.6版本.

① 水力分析模块   HEC-RAS模型可对恒定流或非恒定流进行模拟, 其中非恒定流计算基于质量守恒定律和能量守恒定律(Hydrologic Engineering Center, 2016), 基本控制方程为圣维南(Saint-Venant)方程组, 具体方程如下：

连续方程：

(1)

动量方程：

(2)

式中, A为河道过水断面面积(m²)；t为计算时间(s)；Q为流量(m³·s^-1)；x为河道计算长度(m)；q为单位长度的侧向流量(m³·s^-1)；V为流速(m·s^-1)；g为重力加速度(m·s^-2)；z为水位(m)；S_f为阻力坡降.

② 水质分析模块   水质分析模块包括温度、营养物质模拟.其中营养物模拟主要包括溶解态有机氮、溶解态亚硝酸盐、溶解态有机磷、溶解氧等在水中的迁移转化过程和相互间的物理和化学作用(Hydrologic Engineering Center, 2016).水质模型主要根据质量守恒定律.其控制方程为一维对流扩散方程, 具体形式如下：

(3)

式中, C为污染物浓度(mg·L^-1)；A为河道过水断面面积(m²)；Q为流量(m³·s^-1)；D为纵向离散系数；K为营养物综合降解系数；S_i为源汇项(kg).

在营养物质模拟中, HEC-RAS模型可以通过调节参数, 对于一些模型尚未提供的成分或是复杂模型运行前的简单测试, 用户可选择使用任意组分工具(Arbitrary Constituents Box)进行自定义设置, 其控制方程为：

(4)

式中, C为模拟物质的浓度, 是一个状态变量(mg·L^-1)；K为速率常数, 分为增长速率常数(正值)和衰减速率常数(负值)两种(d^-1).

2.2.2 河道和排污口概化

研究区域水系较为复杂, 模型构建前须对河道进行概化.河道概化要确保概化后的河道输水能力和调蓄能力基本与实际情况相同, 具体操作如下：①骨干河道按原长度进行输入；②汇入骨干河道的相近次要输水河道进行适当合并；③不起输水作用的较小河道适当忽略；④小型湖泊、池塘等概化为零维的调蓄节点.通过以上操作, 最终将研究区域划分为13个河段和93个河段, 保证研究区域主河道两相连断面的间距在300~500 m, 次要河流两相连断面间距不超过1000 m.

同时, 为简化计算, 暂不考虑由于大气沉降、底泥释放以及蒸发沉降等因素对模型的影响, 只考虑点源污染(工业点源、污水处理厂尾水、规模化畜禽养殖)和非点源污染(畜禽分散养殖、水产养殖、农田流失、生活污染等).根据城区地下管网、排污口设置, 结合实际情况将排污口进行概化(图 2).HEC-RAS模型提供了一种将质量而不是浓度引入模型的方法, 通过添加质量注入(Mass Injection)的断面位置、选择物质种类并输入开始时间、持续时间和注入质量, 将模拟组分以质量的形式引入下游的水质单元(Hydrologic Engineering Center, 2016).

2.2.3 模型参数率定

模型参数率定主要针对水动力模拟及水质模拟的相关主要参数进行模拟.其中水动力模拟中的关键参数曼宁系数与水质模拟中的离散系数采用经验估算模拟, 水质衰减系数由于取值范围跨度较大, 采用Matlab进行联合仿真求解.

水动力模拟所需的曼宁系数n先根据实地调研和用户手册推荐值(表 1)进行取值, 再根据结果进行微调, 各河段曼宁系数的最终取值在0.2~0.35之间.

由于目前模型没有专门针对氨氮和总磷的模拟模块, 故采用任意组分的形式进行建模, 所需的主要参数为离散系数(Dispersion Coefficients)和衰减系数(Decay Rate).离散系数D用来表示河流纵向混合情况, 受到河道断面形态、水文条件、河流状态等因素的影响.目前离散系数多利用经验估算, 根据用户手册提供的费舍尔等式(Fischer′s equation)进行计算, 具体形式如下：

(5)

其中，

式中, m为用户自定义乘数, 取1；u为断面平均流速(m·s^-1)；w为河道平均宽度(m)；Y为河道平均深度(m)；u*为剪切速度(m·s^-1)；g为重力加速度, 取9.81 m·s^-2；d为河道平均深度(m)；S为摩擦坡降；n为曼宁系数；R为水力半径, 即过水断面面积与湿周之比(m).

衰减系数k随所在地区、河道情况等改变, 取值范围较大.为减少工作量, 提高参数精确度, 选用Matlab进行联合仿真, 通过循环结构、利用代价函数选择参数最优取值, 具体流程如图 3所示.

按目的主要分为打开和运行程序、更改参数、读取结果、比较和输出4个功能模块, 主要操作如下：

① 打开和运行程序：主要涉及的函数有actxserver函数(创建HEC-RAS软件对应的COM服务器)、h.Project_Open函数(打开工程文件)、h.Compute_CurrentPlan函数(控制程序, 使其根据预设的模拟计划运行)、h.Compute_Complete函数(用于判断运算是否完成).

② 更改参数：先基于copyfile函数, 对水质输出文件*.w##进行镜像；借助fgetl函数逐行读取文件内容, 并使用regexp函数将内容进行分割；使用strcmp函数将内容与目标字符串比较判断；使用fprintf函数将目标数值进行替换.

③ 读取结果：HEC-RAS软件计算结果以二进制的形式存储于*.wq##文件中, 需要借助matlab程序进行解码：使用fseek函数, 令文件指针偏移到数据存储位置, 再借助fread函数, 以“float32”的形式对数据进行解读.

④ 比较和输出：通过之前已整理输入的实测数据(.csv文件), 以平均相对误差为判断依据, 输出相对理想的参数取值.考虑到可能会存在异常值, 利用sort函数进行排序, 去掉5%的异常点再进行比较.其他涉及的函数和语句还有sum函数、min函数和fprintf语句.

2.2.4 模拟效果评估

模型的率定效果及验证效果主要通过纳什系数E_NS、决定性系数R²和相对误差PBIAS来衡量.纳什系数用于衡量实测数据和模拟数据的拟合程度, 决定性系数表示实测数据和模拟数据趋势的一致性, 相对误差表示模拟结果较实测数据的偏移程度.一般来说, 当E_NS≥0.5、R²≥0.6且|PBIAS|≤15%时, 模型模拟结果相对较好(Moriasi et al., 2015).纳什系数E_NS、决定性系数R²和相对误差PBIAS的具体计算公式如下：

(6)

(7)

(8)

式中, O_i为第i个观测值；S_i为第i个模拟值；O_avg为全部观测值的平均值；S_avg为全部模拟值的平均值.

2.3 数据来源

本次研究中污染负荷估算所涉及的人口、畜禽养殖、耕地、水产养殖等数据主要来自于《长兴县统计年鉴》(2017年)；工业点源、城镇生活污染负荷、规模化畜禽养殖污染及污水处理厂等信息主要来自于《长兴县环境统计年鉴》(2016年)；研究区域河道长度、河宽、水深等信息主要来自于《长兴县水域调查报告》(2018年)；2011—2017年河道流量数据主要由浙江省水文局提供；2011—2017年相关断面水质数据主要来自于长兴县环境监测站.

2019年5月对研究区域主要河道、交叉点附近共设置15个监测点进行水文水质同步监测.每天取样3次, 连续监测3 d, 其中流量使用手持式超声波多普勒流速仪和旋杯式流速仪测得, 氨氮、总磷分别根据相关规范、采用纳氏试剂分光光度法、钼酸铵分光光度法进行测试.

3 结果与讨论(Results and discussion) 3.1 水动力参数率定及验证

在模型非恒定流模拟中, 选择流量过程线作为边界条件, 采用2011—2016年水文数据进行率定, 2017年水文数据进行验证.

根据河道情况选取用户手册中的推荐值, 再根据结果进行微调, 最终各个河段的参数取值分布在0.020~0.045之间.水动力率定、验证结果满足模拟精度要求, 其中率定期E_NS=0.99、R²=0.98、PBIAS=2.73%, 验证期E_NS=0.96、R²=0.98、PBIAS=2.9%, 模拟趋势与实际情况一致, 峰值略有偏差(图 4), 整体模拟结果较为理想.

3.2 水质参数率定及验证 3.2.1 水质参数率定

在水质模拟中, 以2012年1—6月为预热期, 以2012年7—2016年6月为率定期.

离散系数根据费舍尔等式计算, 率定选择模型自带的自适应功能(Computed Values), 这一功能已被证实具有较好的适应性和准确性(Szomorova et al., 2015).采用上下浮动15%的范围进行约束, 使参数更好的适应模型, 最终得到各河段的离散系数取值在0.3~1.5之间, 主要河道的离散系数在1~1.5之间, 小于常规河道, 这可能与河道较窄、流速较慢有关.

利用Matlab软件进行联合仿真, 以0.001为步长、在0~0.5的范围内(张亚丽等, 2015)进行衰减系数k的循环尝试, 以提前输入的实测数据为判断依据, 利用代价函数进行计算、排序, 在全部程序运行结束后, 输出最佳的参数取值, 其中k_氨氮=-0.069, k_总磷=-0.044, 这与相近水系的参数率定结果(庄巍等, 2016)相近.率定结果如图 5所示, 模型模拟趋势与实际情况相似, 但峰值模拟效果较差.

对模拟值与实测值进行纳什系数E_NS、决定性系数R²和相对误差PBIAS分析(表 2), 其中氨氮模拟R²＞0.65, E_NS＞0.64, |PBIAS|＜7%；总磷模拟R²＞0.68, E_NS＞0.58, |PBIAS|＜7%, 基本满足E_NS≥0.5、R²≥0.6且|PBIAS|≤15%的要求.

3.2.2 水质模拟验证

以2016年7月—2017年12月3个断面的监测数据(逐月数据)和2019年5月15个点位的实测数据分别进行验证.

2016—2017年3个常规监测断面的水质模拟结果显示, 整体模拟结果的趋势与实测值趋势基本趋于一致(图 6), 但仍存在一定的偏差.根据模拟值与实测值的纳什系数E_NS、决定性系数R²和相对误差PBIAS等分析, 氨氮模拟R²＞0.6, E_NS＞0.53, |PBIAS|＜13%；总磷模拟验证期R²＞0.56, E_NS＞0.52, |PBIAS|＜8%, 基本满足要求, 但相对来说精度还不是很高.

利用2019年5月20—22日的实测数据, 进行多点位、短步长的水质模拟验证, 纳什系数E_NS、决定性系数R²和相对误差PBIAS等计算结果显示, 氨氮模拟R²=0.83, E_NS=0.86, PBIAS=-1.4%；总磷模拟R²=0.81, E_NS=0.85, PBIAS=6.1%, 结果较为理想, 确保了模型在整个研究区域的适应性, 进一步证明了模型的可用性.如图 7也表明了实测值与模拟值趋同性较好.

3.3 误差分析

模型在水动力模拟中结果较好, 趋势和峰值模拟都较为理想；在水质模拟中, 结果基本在可接受范围内, 模拟趋势与实际情况相似, 但峰值存在较大误差, 从两次验证结果来看, 针对短时段的模拟效果远好于长时段的模拟效果.综合分析误差来源的, 主要有以下几方面：

① 模拟时间跨度和模拟步长的选择.时间跨度是影响模型模拟效果的主要因素之一.水质模拟的时间范围是2012—2017年.在这6年时间中, 污染物的排放量、河道的结构和调蓄能力势必会发生改变, 季节环境的变化、政策的实施以及突发事件都会对水质带来影响, 因而在长时段内使用相同的参数(包括污染源概化、污染负荷), 可能会带来较大的误差.第二次验证中, 模拟步长间隔短、模拟时间跨度小, 避免了污染负荷变化、水力调控变化等因素产生的额外影响, 模拟结果相对较好.

② 研究区域地形水文的复杂性.研究区域处于平原地区, 地面起伏不大, 但河网纵横, 又临近太湖, 存在复杂的水量交换问题.入湖口受到太湖的影响, 水流流向不定, 滞流、倒流现象时常发生, 加大了水质模拟的难度.水质率定和第一次验证所选取的断面都是在入湖口附近, 尽管在模型模拟中设置了感潮边界(Tidal Boundary), 仍然不能够较好的降低误差；第二次验证涉及到了研究区域的各个点位, 并且取样期间常有降水, 减少了由于逆流带来的影响, 模拟结果相对较好.

③ 模型模块的不成熟性.现有模型对于水动力以及溶解氧、盐分的模拟相对成熟, 但是还未推出特别针对氨氮和总磷的模拟模块, 模型模拟是以任意组分的方式进行, 只考虑物质在水中的扩散降解作用, 忽略了物质各类的转化作用.

④ 模型概化带来的信息缺失.模型在进行排污口概化、水系概化、污染量核算、边界条件输入时, 很难做到完全与实际情况相同, 同时模型也未考虑底泥与水体交换、大气沉降等因素, 这些都会影响模型模拟结果的精度.

4 结论(Conclusions)

1) HEC-RAS模型常被用于解决水力模拟等相关问题, 近几年新增的水质模拟功能尚未得到广泛应用.本研究结果显示, HEC-RAS模型可以解决交叉/环状河网的模拟问题.同时软件具有体积小、运算速度快, 界面友好, 完全开源等优点, 可以自由下载以及与其他模型/软件联合使用.

2) 模型应用过程中需要对其内含参数进行参数的率定.本研究利用Matlab进行自动化求取参数, 相对于常规的试错法参数率定, 节约了率定时间及人工选择参数的不确定性, 提高了模型参数的率定效率, 也证明了其可兼容性.

3) HEC-RAS在研究区域的水质模拟结果显示, 该模型对于区域环境信息充分、短时段的水质模拟具有较好的精度.然而, 该模型的水质模块仍处于开发阶段, 部分水质成分尚未加入, 配套程序也不是十分齐全, 仍需进一步发展和完善.