河流动态水环境容量核算与影响因素分析

崔 嵩 ，贾朝阳 ，付 强 ，胡 鹏，宋梓菡，卜鑫宇

（1.东北农业大学水利与土木工程学院，哈尔滨 150030；2.东北农业大学松花江流域生态环境保护研究中心，哈尔滨 150030；3.中国水利水电科学研究院，北京 100038）

水环境容量（Water environmental capacity，WEC）是指水体在规定环境目标下所能容纳的最大污染物量，反映水环境在满足特定功能条件下对污染物的承受能力[1-2]。水环境容量是水污染控制和环境管理重要组成，受水系特征、水文条件、污染物理化性质及水质目标等因素影响[3]。目前，在控制污染物排放总量背景下，水环境容量受科研人员关注[4]。20世纪60年代，日本学者首次提出污染物排放总量控制问题后，科研人员针对水环境容量控制问题开展大量研究工作，先后开发QUASAR[5]、MIKE[6]、SWAT[7]、QUAL2E-UNCAS[8]

等模型。自20世纪70年代，我国开始水环境容量研究工作以来，先后建立CJK3D[9]、水动力-水质-水生态耦合模型[10]和水-生态-底泥耦合模型[11]等多种数学模型，目前国内学者在水环境容量核算已取得丰富研究成果[12]。在水环境容量核算过程中，通常应用Monte Carlo 方法研究水环境容量模型不确定性问题[13]。Monte Carlo方法可将水环境容量模型不确定性问题转化为对模型参数的估计量，通过建立不确定性和估计概率间联系，以后验分布概率代表不确定事件发生概率，提升水质管理有效性。

松花江作为我国七大水系之一，是我国“十三五”环境保护规划中重点监管流域，连通哈尔滨、佳木斯、吉林和松原等主要城市，是黑龙江省与吉林省生活和工业用水主要来源，同时也是生活污水和工业废水受纳水体。在哈尔滨段水环境容量核算方面，姜曼曼等采用经济效用模型并结合层次分析法核算化学需氧量（COD）水环境容量[14]。王思文应用WASP模型对松花江哈尔滨段水质现状开展模拟研究，采用试错法计算水环境容量，表明水环境容量遵循丰水期＞平水期＞枯水期的一般规律[15]。马欢通过一维水质模型及二维水质模型对松花江哈尔滨段水环境容量作对比计算，并核算需削减的排污量[16]。

目前，针对水环境容量研究集中于水环境容量核算和分配等方面，而影响核算结果模型中降解系数动态计算及识别模型主要影响参数的不确定性分析方面研究较少。因此，本研究以松花江哈尔滨段为例，以研究区域水文气象数据为驱动因子，综合模型与统计方法，建立降解系数动态计算公式，结合水功能区水质要求与汇入支流水质情况，采用基于二维水质模型的分段求和模型核算松花江哈尔滨段COD、总氮（TN）和总磷（TP）动态水环境容量；同时，通过分析水环境容量计算模型不确定性与输入参数敏感性，揭示影响水环境容量主控环境因子，为人类活动影响下河流环境综合治理及流域经济高质量协同发展提供科学依据。

1 材料与方法

1.1 研究区域概况与控制单元划分

以松花江哈尔滨段（拉林河断面至大顶子山断面，全长约153.5 km）为研究区域，沿途包含运粮河、何家沟、马家沟、阿什河和呼兰河等一级汇入支流。研究区域地处中温带大陆性季风气候区，年平均气温5.6 ℃[17]，水温和流量具有明显季节性特征，每年3 至6 月为平水期，7 至10 月为丰水期，11至次年2月为枯水期（冰封期）[18]。哈尔滨市行政区域总面积53 076.5 km2，建成区面积445.8 km2，建成区外以耕地和林地为主，其中耕地面积约22 706.3 km2，占行政区域总面积42.78%[19]。本研究综合考虑松花江哈尔滨段主要饮用水取水口分布及水体功能要求等综合情况划分控制单元（见图1），其水文特征及水体功能区分布情况见表1。

表1 控制单元划分情况Table 1 Division of control units

图1 研究区域概况Fig.1 Overview of the study area

1.2 数据来源

本研究松花江哈尔滨段水质数据来源于课题组前期研究成果和相关文献[20-21]；水文数据（水温、流量）来源于《中华人民共和国水文年鉴——黑龙江省流域水文资料第3 册》[22]；气象数据（气温、降水量）来源于美国国家海洋和大气管理局（https://www.ncei.noaa.gov/）[23]。

1.3 水环境容量模型与参数选取

1.3.1 水环境容量模型

根据《水域纳污能力计算规程》（GB/T25173-2010）规定[24]，年平均流量大于150 m3·s-1河段属于大型河段，应选取二维水质模型计算纳污能力。松花江哈尔滨段多年平均年径流量为402.3×108m3，折算成年平均流量为1 276 m3·s-1＞150 m3·s-1[25]，因此选取二维水质模型计算其纳污能力。该模型假设污染物排入河流后垂向（z轴方向）瞬间混合均匀，横向y方向和纵向x方向有相对持久浓度梯度。鉴于本研究关注的污染物主要为岸边排放，故采用二维岸边排放模型。

按照上述模型，根据各月份水文与气象条件计算各时段水环境容量。松花江作为哈尔滨市饮用水水源[26]和重要受纳水体，为充分利用松花江哈尔滨段水环境容量，本文采用分段求和模型（Sub⁃section summation model，SSM）核算松花江哈尔滨段水环境容量。SSM根据排水口、取水口和支流等增加控制断面，将功能区划分为若干区段，分别计算各区段水环境容量，再相加得到功能区水环境容量。研究区域控制断面概化如图2所示。

图2 研究区域控制断面概化示意图Fig.2 Schematic diagram of the control section of the study area

假定控制断面前后污染物浓度满足质量守恒定律，将控制单元内支流入口设置为n个，将水功能区划分为n+1个控制区段。根据连续性方程确定每个控制区段流量和污染物浓度，将其作为紧接控制断面后管段初始条件，通过二维水质模型计算下一控制断面前污染物浓度，每个控制区段WEC 之和作为水功能区整体WEC。基于二维水质模型、物质平衡方程和流量平衡方程建立SSM 模型，其方程为：

式中，c（x,y）为污染物在河流中任意敏感点（x,y）的预测质量浓度（mg·L-1）；x为敏感点到排污口纵向距离（m）；y为敏感点到排污口所在岸边横向距离（m）；Qp为排污口流量（m3·s-1）；Cp为排污口污染物排放质量浓度（mg·L-1）；Ch为预测污染物在河流中背景质量浓度（mg·L-1）；K为污染物综合降解系数（d-1）；B为河流水面宽度（m）；H为河流平均水深（m）。

敏感点位置通常取在岸边，故敏感点到排污口所在岸边横向距离为0。上式简化为：

假设各水功能区水质目标浓度为水功能区末端断面水质控制标准，各段水环境容量计算公式为：

式中，Cs为第i个控制区段水质目标浓度（mg·L-1）；xi为第i个控制区段长度（m）；C0为第i个控制区段初始污染物浓度（mg·L-1）。

各功能区段水环境容量之和为：

式中，n为断面总数。

《地表水环境质量标准》（GB3838-2002）给出某一水质等级下水质目标浓度Cs为数值范围[27]，河段水质目标浓度设定无参考标准，导致水环境容量不确定性。对于由n+1个区段组成的水功能区，水环境容量是每个区段水环境容量总和。水环境容量计算将转化为一个带约束的多步决策过程问题，可用优化算法求解。

根据水环境容量定义，水环境容量目标是最大化水环境容量，因此水环境容量数学模型为：

1.3.2 选取水文参数

本文水文参数取值参照姜曼曼等[14]研究，取相应月份最小流量为当月设计流量，其他水文参数与此对应。纵向弥散系数My是反应河流纵向混合输送特性重要参数，其数值与河流水力条件密切相关，受河床粗糙度、河道弯度和河道不规则断面形状等因素影响。本文纵向弥散系数My参照杜慧玲等[28]研究成果。

1.3.3 模拟河流水温

前人研究大多将污染物综合降解系数设为固定值或按不同时期、不同水文条件取值，导致模型模拟结果与真实情况存在较大偏差。在水环境容量计算中，污染物综合降解系数是关键参数，受温度影响较大。本文根据松花江哈尔滨段4～11月历史气温和水温作分析（见图3），提出基于气温的水温经验公式：式中，T为水温（℃）；t为气温（℃）。

图3 气温与水温拟合曲线Fig.3 Fitting curve of air temperature and water temperature

1.3.4 污染物综合降解系数K选取

河流污染物降解过程受水文水力学条件影响较大，因此本文将室内试验法与经验公式法结合计算污染物综合降解系数K。通过现场采集水样，并在东北农业大学国际持久性有毒物质联合研究中心实验室（20 ℃、静态水体）测定各污染物综合降解系数，通过流速、河流坡度、水深和水温等相关因素校正[29]，公式如下：

式中，K20为实验室条件下污染物降解系数（d-1）；u为河流流速（m·s-1）；H为河流水深（m）；α为河床活度系数，无量纲。

由于冰封期无法应用经验公式预测水温，本研究参考李芳圆等[30]研究，12月至次年3月冰下水温依次按0.74 ℃（12月）、0.75 ℃（1月）、0.97 ℃（2月）和1.09 ℃（3月）计算。

2 结果与分析

2.1 理想水环境容量

理想水环境容量是指在忽略支流、排污口和上下游水体影响后，将拉林河-朱顺屯江段（R1）、朱顺屯-东江桥江段（R2）和东江桥-大顶子山江段（R3）视为独立江段，采用段尾控制法计算水环境容量。松花江哈尔滨段3 个控制单元内COD、TN和TP 理想动态水环境容量如图4 所示，COD、TN和TP理想水环境容量最大值出现在7或8月，最小值出现在1、2、6 月。全年水环境容量最大是R3，其COD、TN 和TP 水环境容量分别为155 937.69 t、7 467.86 t、1 514.54 t；COD 环境容量最小为R2，4 309.95 t，TN 和TP 水环境容量最小为R1，分别为3 052.19 t 和643.77 t。TN 与TP 理想水环境容量变化曲线趋势相似，而COD 理想水环境容量则与TN 和TP 存在差异，主要原因为COD 与TN、TP 相比水体容许浓度大，降解系数高，故水环境容量较高。

图4 松花江哈尔滨段理想动态水环境容量Fig.4 Ideal dynamic water environmental capacity in Harbin section of Songhua River

R1的COD、TN 和TP 理想动态水环境容量最大值出现在7月，分别为19 396.41、538.83和114.50 t；COD 的最小值出现在1 月，TN 和TP 的最小值均出现在2 月，属于松花江的冰封期，其理想动态水环境容量分别为每月3 717.60、103.59 和22.05 t；最大值与最小值之比分别为5.22、5.20 和5.19，说明该江段TN、TP 波动幅度小于COD 波动幅度。

R2 的COD 理想动态水环境容量变化曲线较平缓，最大值在8月，为5 073.61 t，最小值在6月，为2 467.68 t，比值为2.06，变化幅度较小，主要受水体流速和水文条件影响，TN与TP理想动态水环境容量最小值均在6月，为189.57和38.42 t，因6月松花江流量较低，水文状况抑制TN和TP降解作用。

R3 的COD 环境容量季节性变化明显，最小值在6月，为8 866.21 t，此时江段流量较低，水体流速较小，水文状况抑制COD降解作用；最大值在8月，为19 184.08 t，此时属松花江哈尔滨段丰水期，降水量较大，水体水文状况利于污染物降解。TN 和TP 环境容量最小值均在6 月，为398.45和82.52 t，最大值均在8 月，为919.70 和186.46 t，最大值与最小值之比分别为2.31和2.26，说明这一江段TN、TP波动幅度较小。

2.2 实际水环境容量

在实际计算水环境容量过程中，需考虑支流及上下游水质对干流水质及水环境容量影响。计算理想水环境容量时通常假设支流水质较好、流量小，忽略该江段支流水质水量对干流影响，但实际监测结果表明支流在水质水量等方面对干流影响不能忽略。因此在计算实际水环境容量时，应在理想水环境容量计算基础上考虑支流汇入情况。计算结果如图5所示。

图5 松花江哈尔滨段实际动态水环境容量Fig.5 Actual dynamic water environmental capacity in Harbin section of Songhua River

实际水环境容量与理想水环境容量在趋势上具有相似性。松花江哈尔滨段COD、TN和TP全年实际水环境容量分别为219 723.26、8 040.78 和2 161.34 t。松花江哈尔滨段实际水环境容量较大。其中，R3实际水环境容量最大，该段COD、TN和TP 实际水环境容量分别占松花江哈尔滨段52%、51%和41%。

松花江哈尔滨段COD、TN 和TP 实际水环境容量最大值在7或8月，最小值在1、2和6月，与理想水环境容量一致。全年实际水环境容量最大为R3，其COD、TN、TP 实际水环境容量为114 206.07、3 320.24 和1 093.66 t；COD 环境容量最小为朱顺屯-东江桥段，为26 303.32 t，TN和TP水环境容量最小为R1，为475.71 和2 143.99 t。丰水期COD 水环境容量最大值（14 325.03 t）、TN 水环境容量最大值（593.66 t）、TP 水环境容量最大值（270.95 t）、平水期COD 水环境容量最大值（10 098.76 t）、TN水环境容量最大值（289.47 t）、TP水环境容量最大值（97.85 t）与枯水期COD 水环境容量最大值（10 401.46 t）、TN 水环境容量最大值（288.84 t）、TP水环境容量最大值（103.54 t）均出现在R3，主要原因是R3江段较长，段首与段尾水质相差较大，容许排放较多污染物。

2.3 季节性变化情况分析

本文以12～2月（冬季）、3～5月（春季）、6～8月（夏季）、9～11 月（秋季）为研究区域季节划分标准，分析2019 年松花江哈尔滨段理想水环境容量与实际水环境容量季节性变化（见图6）。COD、TN和TP 理想水环境容量与实际水环境容量最大值均在夏季或秋季，最小值在冬季和春季，表明松花江哈尔滨段在夏、秋两季对污染物纳污能力较强，而冬季和春季纳污能力较弱。夏季比冬季温度更高，导致夏季活性污染物降解更多，且夏季剧烈降雨引起强烈水交换活动，均可导致水环境容量增加。

图6 松花江哈尔滨段COD、TN、TP的理想与实际水环境容量季节性变化Fig.6 Ideal and actual WEC seasonal variation of COD,TN,and TP in Harbin section of Songhua River

2.4 影响因素分析

河流动态水环境容量主要影响因素源于模型各输入参数（水深H、流量Q、纵向弥散系数My、宽度B和气温T）。本研究对数据概率分布作拟合，基于参数概率密度分布，选用Monte Carlo 模拟方法识别COD、TN和TP水环境容量模型中主要影响因素[31]。应用Crystal Ball 软件对COD、TN 和TP水环境容量作500 000次模拟抽样，敏感性相关等级系数结果见图7，松花江哈尔滨段COD、TN、TP 实际水环境容量主要受气温T影响，其次是流量Q影响，而水面宽度B、水深H和横向扩散系数My对模型影响较小。气温是影响松花江哈尔滨段水环境容量影响主因，是模型反映研究区水环境容量实际情况关键参数。河流水体温度直接影响污染物自净过程中生化作用，水温升高增加生化反应速率和微生物活性，有利于污染物降解，因此在相同条件时水温较高河流综合污染物降解系数高于水温较低河流。

图7 松花江哈尔滨段水环境容量各计算参数敏感度Fig.7 Sensitivity of calculated parameters of water environmental capacity in Harbin section of Songhua River

哈尔滨市全年降水集中在5至9月，约占全年总降水量90.30%。COD、TN和TP实际水环境容量与降水量对应关系如图8 所示，同时Pearson 相关性分析表明，COD 和TP 水环境容量与降水量具有相关性（COD：r=0.766，P=0.004；TP：r=0.63，P=0.028）。长历时水量积累有助于水位上涨，增大河流水量，提高污染物水环境容量。降水量增加，可增加水体稀释能力，促进水体流动和交换，增强水体自净能力，增加水体对污染物容纳能力，因此降水量差异性对水环境容量具有一定影响。

图8 水环境容量与降水量对应关系Fig.8 Corresponding relationship between water environmental capacity and precipitation

3 讨论

水环境容量核算在污染物总量控制、流域水质目标管理及水功能区限制纳污红线管理等方面起重要作用[1]。本研究建立基于气温的水温经验公式，改进综合污染物降解系数计算，在水温数据缺失情况下可求解任意时刻污染物综合降解系数，降低以往因污染物综合降解系数取固定值或按不同时期、不同水文条件取值产生的误差。采用二维水质模型分段求和模型（SSM）考虑支流汇入对水环境容量影响，并通过实际汇入位置计算污染物浓度情况，避免传统计算方法产生误差，可计算水环境容量最大值，保证在容许浓度内水质处于达标状态方法。

本研究以松花江哈尔滨段为例，采用SSM 核算COD、TN和TP水环境容量，其结果小于姜曼曼等使用段尾控制法核算结果[14]。在变化趋势方面，采用SSM核算的水环境容量与杜慧玲[28]、瞿一清等针对松花江及其他类似河流研究结果相似[8]，均表现为丰水期水环境容量最大，平水期次之，枯水期最小，且实际水环境容量与理想水环境容量变化幅度相似，这是因丰水期河流流量大，水温较高，水体水文状况利于污染物降解，同时丰水期降水量较大，水体流动与交换频繁，增强水体自净能力，增大水环境容量。SSM计算结果主要受气温T、流量Q 等影响，与沈晔娜等研究结果类似[32]。因此，研究对水环境容量核算作不确定性分析与计算参数敏感性分析，揭示影响河流动态水环境容量模型中主要影响因素，可根据主要影响因素对污染物排放作调控与优化，保障河流水质满足国家管控标准。

4 结论

a.松花江哈尔滨段COD、TN和TP理想水环境容量是实际水环境容量1.40、1.80和1.36倍，夏秋两季节水环境容量较高，春冬两季节水环境容量较低。气温和流量是影响水环境容量计算主要因素。COD 和TP 实际水环境容量与降水量呈正相关，水环境容量在降水丰富时期有明显提升。

b.SSM通过优化算法计算水环境容量，解决传统模型中功能区内存在多个支流或排污口概化问题，增加水环境容量计算准确性。改进的综合污染物降解系数计算公式可在水温数据缺失地区计算污染物综合降解系数，降低因取值产生的误差。

c.通过调整综合污染物降解系数及相关水文参数，可将SSM 应用于其他类似河流，为其他类似河流水环境容量核算提供思路与方法。与传统模型相比，SSM可在求解过程中，计算任意控制区段水环境容量，为河流环境综合治理提供参考。