河流動態(tài)水環(huán)境容量核算與影響因素分析

崔 嵩 ，賈朝陽 ，付 強 ，胡 鵬，宋梓菡，卜鑫宇

（1.東北農(nóng)業(yè)大學水利與土木工程學院，哈爾濱 150030；2.東北農(nóng)業(yè)大學松花江流域生態(tài)環(huán)境保護研究中心，哈爾濱 150030；3.中國水利水電科學研究院，北京 100038）

水環(huán)境容量（Water environmental capacity，WEC）是指水體在規(guī)定環(huán)境目標下所能容納的最大污染物量，反映水環(huán)境在滿足特定功能條件下對污染物的承受能力[1-2]。水環(huán)境容量是水污染控制和環(huán)境管理重要組成，受水系特征、水文條件、污染物理化性質(zhì)及水質(zhì)目標等因素影響[3]。目前，在控制污染物排放總量背景下，水環(huán)境容量受科研人員關(guān)注[4]。20世紀60年代，日本學者首次提出污染物排放總量控制問題后，科研人員針對水環(huán)境容量控制問題開展大量研究工作，先后開發(fā)QUASAR[5]、MIKE[6]、SWAT[7]、QUAL2E-UNCAS[8]

等模型。自20世紀70年代，我國開始水環(huán)境容量研究工作以來，先后建立CJK3D[9]、水動力-水質(zhì)-水生態(tài)耦合模型[10]和水-生態(tài)-底泥耦合模型[11]等多種數(shù)學模型，目前國內(nèi)學者在水環(huán)境容量核算已取得豐富研究成果[12]。在水環(huán)境容量核算過程中，通常應用Monte Carlo 方法研究水環(huán)境容量模型不確定性問題[13]。Monte Carlo方法可將水環(huán)境容量模型不確定性問題轉(zhuǎn)化為對模型參數(shù)的估計量，通過建立不確定性和估計概率間聯(lián)系，以后驗分布概率代表不確定事件發(fā)生概率，提升水質(zhì)管理有效性。

松花江作為我國七大水系之一，是我國“十三五”環(huán)境保護規(guī)劃中重點監(jiān)管流域，連通哈爾濱、佳木斯、吉林和松原等主要城市，是黑龍江省與吉林省生活和工業(yè)用水主要來源，同時也是生活污水和工業(yè)廢水受納水體。在哈爾濱段水環(huán)境容量核算方面，姜曼曼等采用經(jīng)濟效用模型并結(jié)合層次分析法核算化學需氧量（COD）水環(huán)境容量[14]。王思文應用WASP模型對松花江哈爾濱段水質(zhì)現(xiàn)狀開展模擬研究，采用試錯法計算水環(huán)境容量，表明水環(huán)境容量遵循豐水期＞平水期＞枯水期的一般規(guī)律[15]。馬歡通過一維水質(zhì)模型及二維水質(zhì)模型對松花江哈爾濱段水環(huán)境容量作對比計算，并核算需削減的排污量[16]。

目前，針對水環(huán)境容量研究集中于水環(huán)境容量核算和分配等方面，而影響核算結(jié)果模型中降解系數(shù)動態(tài)計算及識別模型主要影響參數(shù)的不確定性分析方面研究較少。因此，本研究以松花江哈爾濱段為例，以研究區(qū)域水文氣象數(shù)據(jù)為驅(qū)動因子，綜合模型與統(tǒng)計方法，建立降解系數(shù)動態(tài)計算公式，結(jié)合水功能區(qū)水質(zhì)要求與匯入支流水質(zhì)情況，采用基于二維水質(zhì)模型的分段求和模型核算松花江哈爾濱段COD、總氮（TN）和總磷（TP）動態(tài)水環(huán)境容量；同時，通過分析水環(huán)境容量計算模型不確定性與輸入?yún)?shù)敏感性，揭示影響水環(huán)境容量主控環(huán)境因子，為人類活動影響下河流環(huán)境綜合治理及流域經(jīng)濟高質(zhì)量協(xié)同發(fā)展提供科學依據(jù)。

1 材料與方法

1.1 研究區(qū)域概況與控制單元劃分

以松花江哈爾濱段（拉林河斷面至大頂子山斷面，全長約153.5 km）為研究區(qū)域，沿途包含運糧河、何家溝、馬家溝、阿什河和呼蘭河等一級匯入支流。研究區(qū)域地處中溫帶大陸性季風氣候區(qū)，年平均氣溫5.6 ℃[17]，水溫和流量具有明顯季節(jié)性特征，每年3 至6 月為平水期，7 至10 月為豐水期，11至次年2月為枯水期（冰封期）[18]。哈爾濱市行政區(qū)域總面積53 076.5 km2，建成區(qū)面積445.8 km2，建成區(qū)外以耕地和林地為主，其中耕地面積約22 706.3 km2，占行政區(qū)域總面積42.78%[19]。本研究綜合考慮松花江哈爾濱段主要飲用水取水口分布及水體功能要求等綜合情況劃分控制單元（見圖1），其水文特征及水體功能區(qū)分布情況見表1。

表1 控制單元劃分情況Table 1 Division of control units

圖1 研究區(qū)域概況Fig.1 Overview of the study area

1.2 數(shù)據(jù)來源

本研究松花江哈爾濱段水質(zhì)數(shù)據(jù)來源于課題組前期研究成果和相關(guān)文獻[20-21]；水文數(shù)據(jù)（水溫、流量）來源于《中華人民共和國水文年鑒——黑龍江省流域水文資料第3 冊》[22]；氣象數(shù)據(jù)（氣溫、降水量）來源于美國國家海洋和大氣管理局（https://www.ncei.noaa.gov/）[23]。

1.3 水環(huán)境容量模型與參數(shù)選取

1.3.1 水環(huán)境容量模型

根據(jù)《水域納污能力計算規(guī)程》（GB/T25173-2010）規(guī)定[24]，年平均流量大于150 m3·s-1河段屬于大型河段，應選取二維水質(zhì)模型計算納污能力。松花江哈爾濱段多年平均年徑流量為402.3×108m3，折算成年平均流量為1 276 m3·s-1＞150 m3·s-1[25]，因此選取二維水質(zhì)模型計算其納污能力。該模型假設(shè)污染物排入河流后垂向（z軸方向）瞬間混合均勻，橫向y方向和縱向x方向有相對持久濃度梯度。鑒于本研究關(guān)注的污染物主要為岸邊排放，故采用二維岸邊排放模型。

按照上述模型，根據(jù)各月份水文與氣象條件計算各時段水環(huán)境容量。松花江作為哈爾濱市飲用水水源[26]和重要受納水體，為充分利用松花江哈爾濱段水環(huán)境容量，本文采用分段求和模型（Sub?section summation model，SSM）核算松花江哈爾濱段水環(huán)境容量。SSM根據(jù)排水口、取水口和支流等增加控制斷面，將功能區(qū)劃分為若干區(qū)段，分別計算各區(qū)段水環(huán)境容量，再相加得到功能區(qū)水環(huán)境容量。研究區(qū)域控制斷面概化如圖2所示。

圖2 研究區(qū)域控制斷面概化示意圖Fig.2 Schematic diagram of the control section of the study area

假定控制斷面前后污染物濃度滿足質(zhì)量守恒定律，將控制單元內(nèi)支流入口設(shè)置為n個，將水功能區(qū)劃分為n+1個控制區(qū)段。根據(jù)連續(xù)性方程確定每個控制區(qū)段流量和污染物濃度，將其作為緊接控制斷面后管段初始條件，通過二維水質(zhì)模型計算下一控制斷面前污染物濃度，每個控制區(qū)段WEC 之和作為水功能區(qū)整體WEC?；诙S水質(zhì)模型、物質(zhì)平衡方程和流量平衡方程建立SSM 模型，其方程為：

式中，c（x,y）為污染物在河流中任意敏感點（x,y）的預測質(zhì)量濃度（mg·L-1）；x為敏感點到排污口縱向距離（m）；y為敏感點到排污口所在岸邊橫向距離（m）；Qp為排污口流量（m3·s-1）；Cp為排污口污染物排放質(zhì)量濃度（mg·L-1）；Ch為預測污染物在河流中背景質(zhì)量濃度（mg·L-1）；K為污染物綜合降解系數(shù)（d-1）；B為河流水面寬度（m）；H為河流平均水深（m）。

敏感點位置通常取在岸邊，故敏感點到排污口所在岸邊橫向距離為0。上式簡化為：

假設(shè)各水功能區(qū)水質(zhì)目標濃度為水功能區(qū)末端斷面水質(zhì)控制標準，各段水環(huán)境容量計算公式為：

式中，Cs為第i個控制區(qū)段水質(zhì)目標濃度（mg·L-1）；xi為第i個控制區(qū)段長度（m）；C0為第i個控制區(qū)段初始污染物濃度（mg·L-1）。

各功能區(qū)段水環(huán)境容量之和為：

式中，n為斷面總數(shù)。

《地表水環(huán)境質(zhì)量標準》（GB3838-2002）給出某一水質(zhì)等級下水質(zhì)目標濃度Cs為數(shù)值范圍[27]，河段水質(zhì)目標濃度設(shè)定無參考標準，導致水環(huán)境容量不確定性。對于由n+1個區(qū)段組成的水功能區(qū)，水環(huán)境容量是每個區(qū)段水環(huán)境容量總和。水環(huán)境容量計算將轉(zhuǎn)化為一個帶約束的多步?jīng)Q策過程問題，可用優(yōu)化算法求解。

根據(jù)水環(huán)境容量定義，水環(huán)境容量目標是最大化水環(huán)境容量，因此水環(huán)境容量數(shù)學模型為：

1.3.2 選取水文參數(shù)

本文水文參數(shù)取值參照姜曼曼等[14]研究，取相應月份最小流量為當月設(shè)計流量，其他水文參數(shù)與此對應?？v向彌散系數(shù)My是反應河流縱向混合輸送特性重要參數(shù)，其數(shù)值與河流水力條件密切相關(guān)，受河床粗糙度、河道彎度和河道不規(guī)則斷面形狀等因素影響。本文縱向彌散系數(shù)My參照杜慧玲等[28]研究成果。

1.3.3 模擬河流水溫

前人研究大多將污染物綜合降解系數(shù)設(shè)為固定值或按不同時期、不同水文條件取值，導致模型模擬結(jié)果與真實情況存在較大偏差。在水環(huán)境容量計算中，污染物綜合降解系數(shù)是關(guān)鍵參數(shù)，受溫度影響較大。本文根據(jù)松花江哈爾濱段4～11月歷史氣溫和水溫作分析（見圖3），提出基于氣溫的水溫經(jīng)驗公式：式中，T為水溫（℃）；t為氣溫（℃）。

圖3 氣溫與水溫擬合曲線Fig.3 Fitting curve of air temperature and water temperature

1.3.4 污染物綜合降解系數(shù)K選取

河流污染物降解過程受水文水力學條件影響較大，因此本文將室內(nèi)試驗法與經(jīng)驗公式法結(jié)合計算污染物綜合降解系數(shù)K。通過現(xiàn)場采集水樣，并在東北農(nóng)業(yè)大學國際持久性有毒物質(zhì)聯(lián)合研究中心實驗室（20 ℃、靜態(tài)水體）測定各污染物綜合降解系數(shù)，通過流速、河流坡度、水深和水溫等相關(guān)因素校正[29]，公式如下：

式中，K20為實驗室條件下污染物降解系數(shù)（d-1）；u為河流流速（m·s-1）；H為河流水深（m）；α為河床活度系數(shù)，無量綱。

由于冰封期無法應用經(jīng)驗公式預測水溫，本研究參考李芳圓等[30]研究，12月至次年3月冰下水溫依次按0.74 ℃（12月）、0.75 ℃（1月）、0.97 ℃（2月）和1.09 ℃（3月）計算。

2 結(jié)果與分析

2.1 理想水環(huán)境容量

理想水環(huán)境容量是指在忽略支流、排污口和上下游水體影響后，將拉林河-朱順屯江段（R1）、朱順屯-東江橋江段（R2）和東江橋-大頂子山江段（R3）視為獨立江段，采用段尾控制法計算水環(huán)境容量。松花江哈爾濱段3 個控制單元內(nèi)COD、TN和TP 理想動態(tài)水環(huán)境容量如圖4 所示，COD、TN和TP理想水環(huán)境容量最大值出現(xiàn)在7或8月，最小值出現(xiàn)在1、2、6 月。全年水環(huán)境容量最大是R3，其COD、TN 和TP 水環(huán)境容量分別為155 937.69 t、7 467.86 t、1 514.54 t；COD 環(huán)境容量最小為R2，4 309.95 t，TN 和TP 水環(huán)境容量最小為R1，分別為3 052.19 t 和643.77 t。TN 與TP 理想水環(huán)境容量變化曲線趨勢相似，而COD 理想水環(huán)境容量則與TN 和TP 存在差異，主要原因為COD 與TN、TP 相比水體容許濃度大，降解系數(shù)高，故水環(huán)境容量較高。

圖4 松花江哈爾濱段理想動態(tài)水環(huán)境容量Fig.4 Ideal dynamic water environmental capacity in Harbin section of Songhua River

R1的COD、TN 和TP 理想動態(tài)水環(huán)境容量最大值出現(xiàn)在7月，分別為19 396.41、538.83和114.50 t；COD 的最小值出現(xiàn)在1 月，TN 和TP 的最小值均出現(xiàn)在2 月，屬于松花江的冰封期，其理想動態(tài)水環(huán)境容量分別為每月3 717.60、103.59 和22.05 t；最大值與最小值之比分別為5.22、5.20 和5.19，說明該江段TN、TP 波動幅度小于COD 波動幅度。

R2 的COD 理想動態(tài)水環(huán)境容量變化曲線較平緩，最大值在8月，為5 073.61 t，最小值在6月，為2 467.68 t，比值為2.06，變化幅度較小，主要受水體流速和水文條件影響，TN與TP理想動態(tài)水環(huán)境容量最小值均在6月，為189.57和38.42 t，因6月松花江流量較低，水文狀況抑制TN和TP降解作用。

R3 的COD 環(huán)境容量季節(jié)性變化明顯，最小值在6月，為8 866.21 t，此時江段流量較低，水體流速較小，水文狀況抑制COD降解作用；最大值在8月，為19 184.08 t，此時屬松花江哈爾濱段豐水期，降水量較大，水體水文狀況利于污染物降解。TN 和TP 環(huán)境容量最小值均在6 月，為398.45和82.52 t，最大值均在8 月，為919.70 和186.46 t，最大值與最小值之比分別為2.31和2.26，說明這一江段TN、TP波動幅度較小。

2.2 實際水環(huán)境容量

在實際計算水環(huán)境容量過程中，需考慮支流及上下游水質(zhì)對干流水質(zhì)及水環(huán)境容量影響。計算理想水環(huán)境容量時通常假設(shè)支流水質(zhì)較好、流量小，忽略該江段支流水質(zhì)水量對干流影響，但實際監(jiān)測結(jié)果表明支流在水質(zhì)水量等方面對干流影響不能忽略。因此在計算實際水環(huán)境容量時，應在理想水環(huán)境容量計算基礎(chǔ)上考慮支流匯入情況。計算結(jié)果如圖5所示。

圖5 松花江哈爾濱段實際動態(tài)水環(huán)境容量Fig.5 Actual dynamic water environmental capacity in Harbin section of Songhua River

實際水環(huán)境容量與理想水環(huán)境容量在趨勢上具有相似性。松花江哈爾濱段COD、TN和TP全年實際水環(huán)境容量分別為219 723.26、8 040.78 和2 161.34 t。松花江哈爾濱段實際水環(huán)境容量較大。其中，R3實際水環(huán)境容量最大，該段COD、TN和TP 實際水環(huán)境容量分別占松花江哈爾濱段52%、51%和41%。

松花江哈爾濱段COD、TN 和TP 實際水環(huán)境容量最大值在7或8月，最小值在1、2和6月，與理想水環(huán)境容量一致。全年實際水環(huán)境容量最大為R3，其COD、TN、TP 實際水環(huán)境容量為114 206.07、3 320.24 和1 093.66 t；COD 環(huán)境容量最小為朱順屯-東江橋段，為26 303.32 t，TN和TP水環(huán)境容量最小為R1，為475.71 和2 143.99 t。豐水期COD 水環(huán)境容量最大值（14 325.03 t）、TN 水環(huán)境容量最大值（593.66 t）、TP 水環(huán)境容量最大值（270.95 t）、平水期COD 水環(huán)境容量最大值（10 098.76 t）、TN水環(huán)境容量最大值（289.47 t）、TP水環(huán)境容量最大值（97.85 t）與枯水期COD 水環(huán)境容量最大值（10 401.46 t）、TN 水環(huán)境容量最大值（288.84 t）、TP水環(huán)境容量最大值（103.54 t）均出現(xiàn)在R3，主要原因是R3江段較長，段首與段尾水質(zhì)相差較大，容許排放較多污染物。

2.3 季節(jié)性變化情況分析

本文以12～2月（冬季）、3～5月（春季）、6～8月（夏季）、9～11 月（秋季）為研究區(qū)域季節(jié)劃分標準，分析2019 年松花江哈爾濱段理想水環(huán)境容量與實際水環(huán)境容量季節(jié)性變化（見圖6）。COD、TN和TP 理想水環(huán)境容量與實際水環(huán)境容量最大值均在夏季或秋季，最小值在冬季和春季，表明松花江哈爾濱段在夏、秋兩季對污染物納污能力較強，而冬季和春季納污能力較弱。夏季比冬季溫度更高，導致夏季活性污染物降解更多，且夏季劇烈降雨引起強烈水交換活動，均可導致水環(huán)境容量增加。

圖6 松花江哈爾濱段COD、TN、TP的理想與實際水環(huán)境容量季節(jié)性變化Fig.6 Ideal and actual WEC seasonal variation of COD,TN,and TP in Harbin section of Songhua River

2.4 影響因素分析

河流動態(tài)水環(huán)境容量主要影響因素源于模型各輸入?yún)?shù)（水深H、流量Q、縱向彌散系數(shù)My、寬度B和氣溫T）。本研究對數(shù)據(jù)概率分布作擬合，基于參數(shù)概率密度分布，選用Monte Carlo 模擬方法識別COD、TN和TP水環(huán)境容量模型中主要影響因素[31]。應用Crystal Ball 軟件對COD、TN 和TP水環(huán)境容量作500 000次模擬抽樣，敏感性相關(guān)等級系數(shù)結(jié)果見圖7，松花江哈爾濱段COD、TN、TP 實際水環(huán)境容量主要受氣溫T影響，其次是流量Q影響，而水面寬度B、水深H和橫向擴散系數(shù)My對模型影響較小。氣溫是影響松花江哈爾濱段水環(huán)境容量影響主因，是模型反映研究區(qū)水環(huán)境容量實際情況關(guān)鍵參數(shù)。河流水體溫度直接影響污染物自凈過程中生化作用，水溫升高增加生化反應速率和微生物活性，有利于污染物降解，因此在相同條件時水溫較高河流綜合污染物降解系數(shù)高于水溫較低河流。

圖7 松花江哈爾濱段水環(huán)境容量各計算參數(shù)敏感度Fig.7 Sensitivity of calculated parameters of water environmental capacity in Harbin section of Songhua River

哈爾濱市全年降水集中在5至9月，約占全年總降水量90.30%。COD、TN和TP實際水環(huán)境容量與降水量對應關(guān)系如圖8 所示，同時Pearson 相關(guān)性分析表明，COD 和TP 水環(huán)境容量與降水量具有相關(guān)性（COD：r=0.766，P=0.004；TP：r=0.63，P=0.028）。長歷時水量積累有助于水位上漲，增大河流水量，提高污染物水環(huán)境容量。降水量增加，可增加水體稀釋能力，促進水體流動和交換，增強水體自凈能力，增加水體對污染物容納能力，因此降水量差異性對水環(huán)境容量具有一定影響。

圖8 水環(huán)境容量與降水量對應關(guān)系Fig.8 Corresponding relationship between water environmental capacity and precipitation

3 討論

水環(huán)境容量核算在污染物總量控制、流域水質(zhì)目標管理及水功能區(qū)限制納污紅線管理等方面起重要作用[1]。本研究建立基于氣溫的水溫經(jīng)驗公式，改進綜合污染物降解系數(shù)計算，在水溫數(shù)據(jù)缺失情況下可求解任意時刻污染物綜合降解系數(shù)，降低以往因污染物綜合降解系數(shù)取固定值或按不同時期、不同水文條件取值產(chǎn)生的誤差。采用二維水質(zhì)模型分段求和模型（SSM）考慮支流匯入對水環(huán)境容量影響，并通過實際匯入位置計算污染物濃度情況，避免傳統(tǒng)計算方法產(chǎn)生誤差，可計算水環(huán)境容量最大值，保證在容許濃度內(nèi)水質(zhì)處于達標狀態(tài)方法。

本研究以松花江哈爾濱段為例，采用SSM 核算COD、TN和TP水環(huán)境容量，其結(jié)果小于姜曼曼等使用段尾控制法核算結(jié)果[14]。在變化趨勢方面，采用SSM核算的水環(huán)境容量與杜慧玲[28]、瞿一清等針對松花江及其他類似河流研究結(jié)果相似[8]，均表現(xiàn)為豐水期水環(huán)境容量最大，平水期次之，枯水期最小，且實際水環(huán)境容量與理想水環(huán)境容量變化幅度相似，這是因豐水期河流流量大，水溫較高，水體水文狀況利于污染物降解，同時豐水期降水量較大，水體流動與交換頻繁，增強水體自凈能力，增大水環(huán)境容量。SSM計算結(jié)果主要受氣溫T、流量Q 等影響，與沈曄娜等研究結(jié)果類似[32]。因此，研究對水環(huán)境容量核算作不確定性分析與計算參數(shù)敏感性分析，揭示影響河流動態(tài)水環(huán)境容量模型中主要影響因素，可根據(jù)主要影響因素對污染物排放作調(diào)控與優(yōu)化，保障河流水質(zhì)滿足國家管控標準。

4 結(jié)論

a.松花江哈爾濱段COD、TN和TP理想水環(huán)境容量是實際水環(huán)境容量1.40、1.80和1.36倍，夏秋兩季節(jié)水環(huán)境容量較高，春冬兩季節(jié)水環(huán)境容量較低。氣溫和流量是影響水環(huán)境容量計算主要因素。COD 和TP 實際水環(huán)境容量與降水量呈正相關(guān)，水環(huán)境容量在降水豐富時期有明顯提升。

b.SSM通過優(yōu)化算法計算水環(huán)境容量，解決傳統(tǒng)模型中功能區(qū)內(nèi)存在多個支流或排污口概化問題，增加水環(huán)境容量計算準確性。改進的綜合污染物降解系數(shù)計算公式可在水溫數(shù)據(jù)缺失地區(qū)計算污染物綜合降解系數(shù)，降低因取值產(chǎn)生的誤差。

c.通過調(diào)整綜合污染物降解系數(shù)及相關(guān)水文參數(shù)，可將SSM 應用于其他類似河流，為其他類似河流水環(huán)境容量核算提供思路與方法。與傳統(tǒng)模型相比，SSM可在求解過程中，計算任意控制區(qū)段水環(huán)境容量，為河流環(huán)境綜合治理提供參考。