基于數字柵格的河網小流域水質數學模型構建與應用

尹海龍，林夷媛，趙東華，石澤敏

（1.同濟大學環(huán)境科學與工程學院，上海200092；2.中交上海航道勘察設計研究院有限公司，上海200120）

我國南方河網地區(qū)河流縱橫，降雨豐富。隨著污水收集率的不斷提高，在城鎮(zhèn)點源污染排放得到有效治理的情況下，村鎮(zhèn)小流域地區(qū)的點源、面源污染排放成為水環(huán)境治理關注的問題。為此，需要建立耦合點源、面源排放與河道水質的數學模型，評估小流域內各類污染排放造成的水質動態(tài)響應，為污染物減排方案的制定提供科學依據［1-4］。

由于河網模擬范圍較大，大多數情況下只能采用數值方法進行模擬［5］。河道水動力水質模型可以全面描述水體中的水動力—水質過程，但是難以準確獲得面源污染的動態(tài)輸入。因此，非常有必要將降雨產匯流模型和河道水動力水質模型耦合，以此建立河流和流域間的水力聯(lián)系。Liu等［6］集成運用HSPF模型和EFDC模型模擬圣路易斯灣河口的水力水質情況，發(fā)現(xiàn)圣路易斯灣河口對密西西比河的水質具有重大影響；Li等［7］采用SWAT和EFDC耦合方法，研究了長潭水庫1952—2010年間58年降水資料的頻率響應與水質模擬分析；Thompson等［8］將水文模型MIKE SHE和水力學模型MIKE11進行耦合，每個時間段都進行數據的交換和迭代計算，通過這種模擬方式反映了坡面與河道洪泛區(qū)之間的相互作用關系；Narasimhan等［9］基于SWAT和WASP建立了水庫水質評價綜合管理模型，結果表明需要削減35%的TN和TP才能明顯減少水庫葉綠素a的濃度；朱瑤［10］等耦合運用SWAT與WASP模型，模擬苕溪流域面源污染負荷與水質響應情況。然而，物理性分布式水文模型還存在水文過程描述的過參數化問題，對空間信息的精度要求也很大程度上限制了模型的整體效果。

為了克服水文過程過參數帶來的不確定問題并提高空間尺度上的模擬精度，本文采用柵格基礎的數字高程模型（Raster-based Digital Elevation Model，DEM）構建方法，自動提取水系和劃分子流域，以DEM柵格單元為基本模擬單元進行產流和產污計算，精細化描述降雨徑流導致的面源污染排放與受納水體間的空間響應關系，且應用方便。在自主研發(fā)的柵格化徑流污染模型的基礎上耦合開放式的河道水動力水質模型。開放式模型的優(yōu)勢在于有助于克服傳統(tǒng)商業(yè)軟件建模面臨的智能通訊和系統(tǒng)集成性能差的問題，可進一步耦合智能算法，為實現(xiàn)參數自動率定、污染物精準溯源等奠定基礎［11-12］。以常州市武進港小流域為例開展實證研究，對河網水質進行日尺度的動態(tài)模擬和精準預測。

1 模型的基本原理

1.1 水動力水質模型

1.1.1 一維非恒定流模型

對于河網地區(qū)而言，河道縱向長度遠大于其寬度和深度，可將河網概化為一維問題。根據質量守恒和動量守恒定律，可建立描述河網水動力過程的一維圣維南方程組，如下所示：

式中：z為河道水位；x為河道縱向距離；t為時間；Q為斷面過流流量；A為橫斷面面積；Bt為河道橫斷面的非流動調蓄面積；ql為單位長度上的旁側入流量，包括點源產生的污水排放量和降雨徑流入河量，其中降雨徑流入河量由柵格化降雨徑流污染模型計算確定，如1.2節(jié)所述；g為重力加速度；R為河道橫斷面的水力半徑；n為河道粗糙系數。

1.1.2 水質模型

河網水質模型控制方程由如下公式描述：

式中：φ為水質組分質量濃度；E為縱向離散系數；SE為點源和面源污染負荷排放匯入，其中面源污染負荷由柵格化面源模型計算確定，如1.2節(jié)所述；SI是描述水體中生物化學反應動力學過程，其作用機理架構如圖1所示。

圖1 水質模擬組分的耦合作用關系描述Fig.1 Description of coupling effect of components in water quality simulation

在圖1中，水質模型描述了水體溶解氧平衡為核心的模型組分之間的相互作用關系。溶解氧的消耗考慮了含碳有機物生化降解耗氧和氨氮硝化作用耗氧以及河道底泥有機質分解耗氧；溶解氧的補充考慮了大氣向水體中的復氧過程。其中氨氮的轉化又涉及到氮的循環(huán)過程，即有機氮、氨氮、硝酸鹽之間的轉化。因此，水質模擬組分主要包括碳化BOD（CBOD）、有機氮、氨氮、硝酸鹽和溶解氧及其之間的降解、轉化速率等。

考慮到水動力和水質模型的開放性以及將來進一步建立智慧水務平臺的需要，水動力和水質模型的模擬計算由開放式的模擬工具HEC-RAS模型系統(tǒng)（Hydrologic Engineering Center’s River Analysis System）實現(xiàn)［13-14］。在此基礎上實現(xiàn)與自主開發(fā)的柵格化降雨徑流污染模型的集成，滿足河網地區(qū)水質動態(tài)模擬和水質動態(tài)達標分析的功能。

1.2 柵格化降雨徑流污染模型

柵格化降雨徑流污染模型以數字等高模型（Digital Elevation Model，DEM）為基礎，自動提取水系和劃分子流域，然后以子流域為基礎計算降雨徑流產流、產污量，并最終分配到對應的河道計算單元。模型的主要功能包括柵格降雨量插值計算提取、柵格徑流深度計算、柵格徑流污染物量計算和柵格水流匯流計算4個模塊，計算思路如圖2a所示。

圖2 柵格化降雨徑流污染模型原理Fig.2 Principles of precipitation-runoff nonpoint pollution model based on DEM raster

（1）降雨量柵格分布計算。由于降雨量站點的實測數據僅為站點所處位置的降雨量，采用反向距離加權法（IDW）進行插值獲得整個流域內每個柵格的降雨量［15］。

（2）降雨徑流深度柵格分布計算。采用SCS曲線值法計算對應不同用地類型的柵格徑流深度。SCS法為降雨徑流關系法，能反映不同土壤類型、不同土地利用方式及前期土壤含水量對降雨徑流的影響，具有參數少、簡單易行等特點［16］。SCS曲線值法的方程為

其中

式中：S為潛在最大入滲量；P為一次降雨總量；Q為實際徑流深；CN為該日的曲線值；Ia為由于地表儲存、截留和下滲而導致的地表徑流的初始損失量，與S呈一定的正比關系，通?？扇ˇ藶?.2。

（3）降雨徑流污染物柵格分布計算。進一步，柵格像元上的降雨徑流攜帶污染物量采用如下公式計算：

式中：L是柵格像元的污染物排放量；Cland為柵格像元的對應土地利用類型污染物徑流質量濃度。

（4）柵格徑流和污染物匯流計算。采用8點重力流向算法（deterministic eight-node，D8），基于數字地形高差進行匯流流向判斷，自動提取水系，劃分匯水區(qū)域［17］。為避免匯流過程模擬的失真，確定水流方向前需要進行DEM數據無凹陷處理。D8算法的基本思想是：在3×3窗口中，計算中心網格與鄰域8個格網之間的距離權落差，具有最大權落差值的鄰域格網方向即定義為水流方向，并且規(guī)定一個柵格單元的水流方向用一個特征碼表示，如圖2b所示。距離權落差計算公式如下：

式中：Δhi是中心柵格與鄰域柵格的高程差值，柵格間的距離d與方向有關，在對角線方向為2倍柵格間距，在其他方向為1倍柵格間距。

根據流向計算結果分析每一柵格單元對應的上游匯水面積值。通過從上游到下游的匯水面積累加，建立河道計算單元與匯水面積的對應關系。結合柵格徑流量和污染物計算結果，將降雨徑流污染分配到對應的河道計算單元。

自主開發(fā)的降雨徑流污染模型優(yōu)勢在于：①產流計算基于柵格式SCS模型，同時考慮了降雨和下墊面條件的空間不均勻性，具有空間尺度上計算精度高的特點；②產污計算基于輸出系數模型，避免了面源污染形成機理的復雜性，所需參數少、操作簡便。

2 模型的建立與率定

2.1 研究區(qū)域概況

武進地處長三角地理中心，位于常州、無錫及宜興三市交界處。武進港是太湖梅梁灣的主要入湖河流之一，是連通京杭大運河和太湖的主要水上通道。武進港全長29 km，水深2～3 m，河寬25～30 m，流入太湖水量枯水年為2.72×108m3·a-1，豐水年為4.36×108m3·a-1。

本文將模型從武進港干流擴展到整個流域范圍內的河網地區(qū)，研究區(qū)域土地利用類型以農業(yè)用地、居住用地、工業(yè)用地為主。概化河道包括武進港及武進港兩岸采菱港、武南河、永安河、錫漂漕河、雅浦港5條主要河流，另外還有禮嘉大河、小留河、虎臣河、政平河、東陽岸河等28條支流。武進港及其流域概況如圖3所示。

圖3 研究區(qū)域概況Fig.3 Description of the site studied

2.2 模型邊界條件

水動力模型上游和下游邊界均采用水位過程線，數據來自坊前、常州、百瀆口等水文站2018年日監(jiān)測數據。

水質模型邊界條件包括點源和面源污染負荷的輸入，點源又可分為集中式點源和分散式點源。研究區(qū)域內集中式點源主要包括7座城鎮(zhèn)污水處理廠尾水排放。分散式點源主要包括未經處理直排的農村生活污水和畜禽養(yǎng)殖廢水，依據《太湖流域主要入湖河流水環(huán)境綜合整治規(guī)劃編制技術規(guī)范》［18］，采用排污系數法對分散式點源的化學需氧量（COD）、氨氮（NH3-N）、總氮（TN）、總磷（TP）負荷量進行估算，計算公式如下：

式中：W為分散式點源污染負荷量；Np為人口數；NL為畜禽養(yǎng)殖量；αp為農村生活污水排污系數，其中COD取27 g·人-1·d-1，NH3-N取4 g·人-1·d-1，TN取6 g·人-1·d-1，TP取0.2 g·人-1·d-1；αL為畜禽養(yǎng)殖排污系數，換算成生豬計算，COD取40.55 g·頭-1·d-1，NH3-N取1.31 g·頭-1·d-1，TN取3.05 g·頭-1·d-1，TP取0.38 g·頭-1·d-1；βP為農村生活污水入河系數，取0.7；βL為畜禽養(yǎng)殖入河系數，取0.6。

面源邊界條件輸入采用降雨徑流污染模型的模擬結果，降雨徑流污染模型中污染排放系數取值如表1所示。需要注意的是，污染負荷統(tǒng)計參數COD需與水質模型輸入參數CBOD進行換算，換算系數依據經驗取值：取5日生化需氧量BOD5與碳的生化需氧量CBOD的比值為0.70，BOD5與COD的比值為0.50，可 推 算 得CBOD與COD的 比 值為0.71［19］。

表1 研究區(qū)域不同土地利用類型徑流面源污染輸出系數Tab.1 Runoff nonpoint pollution export coefficients from different land use types in the catchments studied

計算得研究區(qū)域排放總流量為1.6×108t·a-1，COD、NH3-N、TN、TP污染負荷排放總量分別為9 024.1 t·a-1、621.4 t·a-1、1 275.8 t·a-1、80.8 t·a-1，圖4給出了集中式點源、分散式點源和徑流面源污染的負荷排放占比。可見雖然分散式點源排放流量很小，但排放污染負荷量占比卻很高，COD、NH3-N、TN、TP占 比 分 別 為42.1%、46.3%、58.5%、45.2%。對于NH3-N，降雨徑流污染占比也較高。

圖4 研究區(qū)域各污染類型組成比例Fig.4 Discharge and proportion of each pollution type in the catchments studied

圖5進一步展示了武進港流域以高精度DEM 5m分辨率柵格為單位的氨氮產生量以及不同河段的污染物空間分布，可以看出，雨季（8月）的降雨徑流污染負荷入河量明顯高于旱季（1月）。

2.3 模型率定驗證

水質模型率定包括以下幾個步驟：①輸入上下游水動力邊界及旁側入流量，通過圣維南方程建立河道水動力模型；②輸入點源和面源污染負荷以及各水質參數以模擬河道中水質的時空變化；③比較河道斷面水質的實測值和模擬值，必要時調整水質參數。

分別采用納什效率系數（Nash-Sutcliffe efficiency，NSE）和決定系數（R2）來評價水動力模型和水質模型質量的好壞，計算公式如下：

式中：Qsimi、Qobsi分別為日均流量模擬值和實測值；為日均流量實測值的平均值；n代表總共匹配的天數。NSE取值為負無窮至1，NSE越接近1，表示模型質量越好，模型可信度高。通常認為NSE達到0.75以上代表模擬結果很理想，NSE值在0.75～0.36之間模擬結果也是可以接受的，當NSE低于0.36代表模擬結果不理想［20］。

式中：fsimi、fobsi分別為月均水質濃度模擬值和月水質濃度實測值分別為月均水質濃度模擬值和實測值的平均值；m代表計算時間段內的月份數。R2取值為0至1，R2越接近1，表示模擬精度越好。通常認為R2大于0.6代表模擬結果理想，R2在0.6～0.5之間認為模擬結果總體上可以接受［21］。

研究區(qū)域河網水動力模型采用黃埝橋水文站2018年1—12月的日監(jiān)測流量數據進行率定，如圖6所示?？梢钥闯?，模擬值與實測值吻合情況較好，計算得NSE值為0.71，模擬誤差在允許范圍之內，率定后的河道粗糙系數取值為0.030～0.035。

圖6 黃埝橋斷面流量率定結果Fig.6 Modeling calibration result for flow discharge at Huangnian Bridge station

水質模型采用武進港月采樣監(jiān)測數據（2018年4—12月），率定溶解氧（DO）、氨氮（NH3-N）、高錳酸鹽（CODMn）和總磷（TP）4個水質指標，同樣水質監(jiān)測指標CODMn需與水質模型中的率定參數CBOD進行換算，取BOD5/CBOD=0.70、BOD5/CODMn=0.89，可推算得：CODMn/CBOD=0.79［22］。選取慈瀆橋（CDQ2）和武南河下游（WN7）2個典型斷面進行年尺度模擬。CDQ2斷面位于武進港上游，主要污染源為上游客水，未受支流匯入污染影響；WN7斷面位于武南河和武進港交匯處，水質受采菱港、武南河等匯入支流污染影響，總體而言水質差于CDQ2斷面。由率定結果（圖7）可以看出，模擬值與實測值吻合情況較好。依據式（9）計算得：①針對CDQ2斷面，DO、NH3-N、CODMn和TP的R2值分別為0.82、0.69、0.61和0.71；②針對WN7斷面，DO、NH3-N、CODMn和TP的R2值 分 別 為0.88、0.63、0.59和0.66，說明模擬誤差均在允許范圍之內，水質模型中主要參數率定結果如表2。率定結果顯示，降雨強度較大時，受降雨徑流匯入負荷的沖擊，河道斷面水質短時間內會出現(xiàn)急劇升高的現(xiàn)象，降雨結束后又恢復至常態(tài)。

圖7 CDQ2、WN7斷面水質參數率定結果Fig.7 Calibration results of parameters for water quality modeling at CDQ2 and WN7 station

表2 水質模型主要參數率定結果Tab.2 Calibrated results of main parameters for water quality modeling

3 水質改善方案評估

3.1 污染削減方案比選

針對流域內分散點源占比高的特點，提出優(yōu)先削減分散實點源污染排放，分析污染削減效果與入湖考核斷面之間的水質響應關系。

由2.2節(jié)討論可知，雖然分散式點源排放流量很小，但排放污染負荷量占比卻最高。因此優(yōu)化考慮農村分散式點源污染的削減方案，并評估污染削減前后的河道水質動態(tài)達標頻次。

3.2 畜禽養(yǎng)殖污染控制

通過全面清理整頓非法和不符合規(guī)范標準的養(yǎng)殖場，進行養(yǎng)殖專業(yè)化，可實現(xiàn)畜禽養(yǎng)殖接近“零排放”。模擬研究片區(qū)內畜禽養(yǎng)殖全部關?；蛲ㄟ^還田再利用實現(xiàn)畜禽污染零排河，即片區(qū)內畜禽養(yǎng)殖污染負荷削減率100%的情況下姚巷橋斷面水質的響應情況，模擬結果見圖8和表3。

圖8 畜禽養(yǎng)殖負荷削減前后入湖斷面水質響應情況Fig.8 Predicted water quality at lake inflow station before and after livestock pollution reduction

削減畜禽養(yǎng)殖負荷后水質得到了一定程度的提升，特別是CODMn改善顯著，年達標頻次（地表水III類標準，簡稱III類）從23.0%提升至64.9%，但是NH3-N仍無法達標。

3.3 農村污水處理

收集未經處理的農村生活污水，經處理后以一級A標準排放。在畜禽養(yǎng)殖“零排放”基礎上，結合流域內污染物空間分布特征，比選出對水質影響大的污染優(yōu)先控制村落群，優(yōu)先削減該村落群80%農村生活污水，以期進一步提高水質年達標頻次。

共設計了3種污染優(yōu)先控制村落比選方案：①工況1。選取片區(qū)內農村生活污水排放量最高的區(qū)域7個村落為污染優(yōu)先控制村落，共削減COD 2 591.7t·a-1、NH3-N 120.4 t·a-1、TP 12.9 t·a-1；②工況2。在工況1的基礎上，增加排入武進港干流的10個 村 落，共 削 減COD 3 330.4t·a-1、NH3-N 178.9 t·a-1、TP 30.3 t·a-1；③工況3。在工況1的基礎上，增加排入上游來水京杭大運河的7個村落，共削 減COD 3 205.4 t·a-1、NH3-N 170.9 t·a-1、TP 29.1 t·a-1。分別模擬3種工況下姚巷橋斷面水質的響應情況，模擬結果見圖9和表3。

圖9 不同農村污水削減工況下入湖斷面水質響應情況Fig.9 Predicted water quality at lake inflow station of different rural sewage reduction schemes

表3 不同負荷削減工況下姚巷橋斷面水質達標情況分析Tab.3 Annual attainment rates of NH3-N,CODMn,and TP at Yaoxiang Bridge station of different load reduction schemes

模擬結果可知，進一步削減污染優(yōu)先控制村落相應農村生活污水后，姚巷橋斷面水質得到明顯改善。針對工況1，CODMn和TP改善效果較明顯，CODMn年均值降至5.86 mg·L-1，年達天數（III類）從84 d提升至252d；TP的年均值降至0.22 mg·L-1，年達天數（III類）也從60 d提升至142 d，但是NH3-N仍無法達到地表水III類標準。工況2、3在工況1基礎上分別削減排入武進港和京杭大運河村落群的入河負荷，針對CODMn和TP，工況2、3模擬結果相差不大；但是針對NH3-N，工況3的模擬結果明顯優(yōu)于工況2，NH3-N年均值降低至1.06 mg·L-1，年達標頻次（III類）達66.3%，年達標天數（III類）達242 d。綜上可得工況3為相對最優(yōu)工況。

綜上，改善客水水質、加強農村生活污水處理設施建設以及分散畜禽養(yǎng)殖污染治理可有效改善武進港水質；其中改善上游來水京杭大運河對提高姚巷橋斷面年達標率效果尤為顯著。相對最優(yōu)工況下仍存在不達標天數主要是因為雨天受降雨徑流匯入負荷的沖擊，河道斷面水質短時間內會出現(xiàn)急劇升高的現(xiàn)象。因此從水質動態(tài)達標的角度考慮，若要實現(xiàn)考核斷面全面、穩(wěn)定達標，需進一步加強對農田徑流和地表徑流污染的控制，如推進面源污染的河道生態(tài)攔截等工程措施，從而確保在雨后一段時間內能盡快恢復水質。

4 結論

建立了適用于河網小流域水環(huán)境系統(tǒng)模擬的水文—水動力—水質耦合數學模型，并在常州市武進港流域進行了實證研究。

（1）建立的耦合模型中，降雨徑流污染模型采用自主開發(fā)的模型系統(tǒng)，以DEM柵格單元為基本模擬單位，可以對降雨徑流導致的面源污染時空分布和匯流去向進行精準模擬。在此基礎上，河道水動力水質模型采用開源模型HEC-RAS，構建了流域污染物排放與河道水質的響應關系，實現(xiàn)河網地區(qū)日尺度的水質模擬和動態(tài)預報。模型系統(tǒng)同時具有開放性，能夠為小流域智慧監(jiān)管平臺提供技術支撐。

（2）針對武進港小流域分散式點源污染負荷高、部分匯入支流水質較差、長期難以達標的問題，將水動力水質模型從武進港干流擴展到整個流域范圍內的河網地區(qū)，從而能夠對武進港流域內的污染總量控制進行定量評估。結合流域內入河污染源組分和空間分布特點，模擬不同負荷削減工況下入湖考核斷面姚巷橋的水質響應情況。從動態(tài)達標角度分析考慮，對不同水質方案進行評估。最優(yōu)工況下，針對水質指標CODMn、NH3-N和TP，考核斷面姚巷橋的年達標頻次（III類）分別從23.0%、0、16.4%升高至71.8%、66.3%、75.9%。

作者貢獻說明：

尹海龍：數學模型建立、論文撰寫。

林夷媛：數學模型建立、論文撰寫。

趙東華：技術和材料支持。

石澤敏：技術和材料支持。