喻 婷，陳曉群，苗 滕

(湖北省水利水電科學(xué)研究院 湖北省水利水電科技推廣中心，武漢 430070)

2002年《中華人民共和國(guó)水法》首次明確了“水域納污能力”的概念，此后水域納污能力與水域限制排污總量一起構(gòu)成了我國(guó)水資源保護(hù)行業(yè)的重要基礎(chǔ)，為相關(guān)行業(yè)水資源管理提供了有效依據(jù)[1]。河流納污能力是判斷河流水質(zhì)狀況，提升水資源管理水平的有效方法，國(guó)內(nèi)學(xué)者們根據(jù)其研究河道的特點(diǎn)選取了不同的水質(zhì)模型并進(jìn)行了相應(yīng)的納污能力計(jì)算研究。陶淑蕓等采用MIKE11軟件建立薔薇河水文水動(dòng)力和水質(zhì)模型，計(jì)算和率定薔薇河流域污染負(fù)荷和水功能區(qū)納污能力[2]，但該研究沒(méi)有將污染源進(jìn)行分區(qū)計(jì)算，不利于污染物入河總量控制和管理。閆峰陵等根據(jù)采用二維水質(zhì)模型計(jì)算金沙江攀枝花河段梯級(jí)建設(shè)前后COD、氨氮納污能力，分析納污能力受梯級(jí)開(kāi)發(fā)影響程度[3]，但該研究沒(méi)有優(yōu)選模型參數(shù)及介紹河道設(shè)計(jì)水文條件，計(jì)算結(jié)果合理性亦無(wú)分析。詹曉群分別采用一維、二維非穩(wěn)態(tài)模型和湖庫(kù)均勻混合衰減模型計(jì)算東江源區(qū)水功能區(qū)納污能力，核定限排總量[4]，但未分析污染物削減的類型，不利于治污工作的開(kāi)展。滕昱采用河流一維水質(zhì)模型計(jì)算大連市重點(diǎn)河流的納污能力及污染物控制排放總量[5]，但該研究沒(méi)有按照水功能區(qū)進(jìn)行計(jì)算，且分析評(píng)價(jià)尺度較粗。

本文針對(duì)舉水干流中下游麻城～柳子港段進(jìn)行計(jì)算分區(qū)，分析得到各分區(qū)入河污染源，建立MIKE二維水動(dòng)力對(duì)流擴(kuò)散模型模擬得到設(shè)計(jì)水文條件下河段內(nèi)COD、NH3-N、TP的變化趨勢(shì)并優(yōu)選出水動(dòng)力水質(zhì)模型參數(shù)，最后采用河流一維模型計(jì)算出舉水干流各計(jì)算分區(qū)的納污能力。通過(guò)對(duì)比分析各分區(qū)納污能力與入河污染物量，研究得到各分區(qū)污染物限排方案，為舉水流域開(kāi)展污染整治行動(dòng)提供數(shù)據(jù)支持。

1 水質(zhì)模型建立

舉水干流中下游麻城～柳子港段長(zhǎng)116.69 km，是接納麻城市、武漢市新洲區(qū)工業(yè)廢水和城市生活污水的主要載體。該河段河床起伏不大，地勢(shì)平緩，多為丘陵和河谷平原。為較好地反映舉水污染物空間分布狀況，選用MIKE21模型對(duì)水質(zhì)狀況進(jìn)行模擬。利用實(shí)測(cè)水文資料對(duì)模型計(jì)算所需基本參數(shù)、水動(dòng)力參數(shù)和對(duì)流擴(kuò)散參數(shù)進(jìn)行率定選擇。

1.1 模型介紹

舉水中下游段二維水動(dòng)力及水質(zhì)的模擬分別通過(guò)HD模塊和AD模型進(jìn)行計(jì)算。HD模塊為二維水動(dòng)力模型，其控制方程為二維淺水方程，包括連續(xù)性方程和動(dòng)量方程。AD模塊為對(duì)流擴(kuò)散模型，與水動(dòng)力模塊是動(dòng)態(tài)連接的。模擬系統(tǒng)是基于Boussinesq和靜壓假定的二維不可壓縮雷諾平均Navier-Stokes方程的數(shù)值算法。

1.1.1 水動(dòng)力模型原理

(1)水動(dòng)力模型基本控制方程。笛卡爾坐標(biāo)系下的二維淺水方程：

(1)

(2)

(3)

該二維淺水方程基于Bousinesq渦黏假定和靜壓假定。

(2)湍流模型。湍流建模采用大渦模擬方法中的Smagorinsky壓格子尺度模型。該模型用一個(gè)與特征長(zhǎng)度尺度相關(guān)的有效渦黏值來(lái)描述亞網(wǎng)格尺度輸移。亞網(wǎng)格尺度渦黏值由下式給出：

(4)

式中：cs為定值；l為特征長(zhǎng)度；Sij為形變率。

形變率Sij由下式給出：

(5)

(6)

(7)

曼寧系數(shù)M可以由底床糙率長(zhǎng)度得到：

(8)

1.1.2 水質(zhì)模型原理

MIKE21對(duì)流擴(kuò)散模型根據(jù)HD模塊產(chǎn)生的水動(dòng)力條件，應(yīng)用對(duì)流擴(kuò)散方程進(jìn)行計(jì)算。

對(duì)流擴(kuò)散方程：

(9)

對(duì)流擴(kuò)散系數(shù)是一個(gè)綜合參數(shù)項(xiàng)，包含了分子擴(kuò)散、湍流擴(kuò)散以及剪切擴(kuò)散效應(yīng)。

1.2 計(jì)算分區(qū)及控制因子

舉水流域計(jì)算分區(qū)綜合參考水功能區(qū)劃成果、縣級(jí)行政區(qū)劃等因素，共分為10個(gè)區(qū)?？紤]到污染源經(jīng)概化處理后，計(jì)算工作量會(huì)大幅度減小[6]，本研究分析計(jì)算出各分區(qū)的概化污染源，見(jiàn)表1。

根據(jù)舉水流域水質(zhì)現(xiàn)狀評(píng)價(jià)結(jié)論，水體表現(xiàn)出以有機(jī)類、營(yíng)養(yǎng)類污染為主的典型特征，與流域內(nèi)居民生產(chǎn)生活方式特征相吻合[7]。本研究選取COD、NH3-N、TP作為納污能力分析計(jì)算的控制因子。

1.3 地形文件

由舉水流域麻城～柳子港段地形圖，提取高程點(diǎn)數(shù)據(jù)，生成50×50 m網(wǎng)格地形文件，見(jiàn)圖1。

表1 計(jì)算分區(qū)及概化污染源統(tǒng)計(jì)

Tab.1 Computation of zoning and generalized pollution source statistics

一級(jí)功能區(qū)名稱二級(jí)功能區(qū)名稱縣級(jí)行政區(qū)計(jì)算分區(qū)編號(hào)廢污水量/(萬(wàn)m3·a-1)主要污染物入河量/(t·a-1)CODTPNH3-N舉水源頭保護(hù)區(qū)麻城市P1291.30992.7314.06142.50舉水麻城上游保留區(qū)麻城市P2453.041378.7319.57201.45舉水麻城開(kāi)發(fā)利用區(qū)舉水麻城飲用水源、工業(yè)用水區(qū)麻城市P3144.93483.946.6680.26舉水麻城排污控制區(qū)麻城市P4156.26456.326.2775.86舉水麻城過(guò)渡區(qū)麻城市P5296.30954.7712.94141.65舉水閔家集～三店保留區(qū)麻城市P6732.573140.0341.38459.14新洲區(qū)P7367.941282.2417.67198.64舉水新洲開(kāi)發(fā)利用區(qū)舉水新洲飲用水源、工業(yè)用水區(qū)新洲區(qū)P8268.94688.649.28100.48舉水新洲～團(tuán)風(fēng)保留區(qū)新洲區(qū)P9617.552330.2428.81310.20團(tuán)風(fēng)縣P10461.541984.1622.47308.63合 計(jì)3790.3813691.80179.112018.81

1.4 計(jì)算流量和下游水位

為滿足設(shè)計(jì)條件，在現(xiàn)有資料的情況下，上斷面流量過(guò)程采用2014-2016年麻城水文站日平均流量，下斷面水位過(guò)程采用同時(shí)期柳子港水文站日平均水位，模擬時(shí)段2014年1月至2016年3月。

1.5 排污口參數(shù)

研究河道流經(jīng)麻城城區(qū)，將麻城城區(qū)點(diǎn)源污染排放概化為一個(gè)源點(diǎn)(圖2)，入河流量為0.81 m3/s，COD、TP和NH3-N的濃度分別為47.71、0.85和8.20 mg/L。

圖2 模型計(jì)算區(qū)域源點(diǎn)示意圖Fig.2 Schematic diagram of source points in model

2 模型參數(shù)率定

MIKE 21FM水動(dòng)力水質(zhì)模型參數(shù)主要包括渦流黏度、糙率系數(shù)、擴(kuò)散系數(shù)和衰減系數(shù)。本次研究參考類似河道水動(dòng)力水質(zhì)模擬參數(shù)取值范圍，通過(guò)不斷調(diào)整參數(shù)取值使水質(zhì)模擬值和實(shí)測(cè)值最接近，從而確定模型參數(shù)。

3 模擬計(jì)算結(jié)果

(1)水動(dòng)力模擬情況。由水動(dòng)力模型計(jì)算可知，研究河段內(nèi)研究期內(nèi)流速最大值為0.73 m/s，最小值0.01 m/s，平均流速0.06 m/s(見(jiàn)圖3)。

圖3 研究河道中下游流速模擬值Fig.3 Simulation value of flow velocity in the middle and lower reaches of the river

(2)水質(zhì)模擬情況。通過(guò)水動(dòng)力水質(zhì)耦合模擬，得到研究期內(nèi)河道COD、TP、NH3-N的平均濃度分別為8.39、0.47、0.89 mg/L。通過(guò)與舉水干流麻城水文站、新洲柳子港水文站監(jiān)測(cè)資料以及2018年5月現(xiàn)場(chǎng)實(shí)測(cè)水質(zhì)結(jié)果對(duì)比，模擬結(jié)果能夠較好研區(qū)域內(nèi)水質(zhì)變化規(guī)律。

在日變化過(guò)程中，6、7月污染物濃度顯著大于河流水質(zhì)管理目標(biāo)(見(jiàn)圖4～圖6)。主要因?yàn)榱饔騼?nèi)農(nóng)業(yè)活動(dòng)較多，每年梅雨季節(jié)流域內(nèi)降水顯著增加會(huì)攜帶農(nóng)業(yè)面源大規(guī)模入河，導(dǎo)致水質(zhì)惡化。

圖4 COD模擬濃度Fig.4 Simulated concentration of COD

圖5 TP模擬濃度Fig.5 Simulated concentration of TP

圖6 NH3-N模擬濃度Fig.6 Simulated concentration of NH3-N

4 納污能力計(jì)算

4.1 計(jì)算公式

在納污能力計(jì)算中，參數(shù)的確定和取值是否符合客觀實(shí)際，直接關(guān)系到計(jì)算結(jié)果是否準(zhǔn)確合理。本研究直接采用MIKE21模型率定得到的水動(dòng)力水質(zhì)參數(shù)(表2)進(jìn)行納污能力核算。

表2 水動(dòng)力水質(zhì)模型參數(shù)

Tab.2 Parameters in the Hydrodynamic water quality model

參數(shù)名稱取值說(shuō)明模擬時(shí)段20141201至20160301舉水中下游段低水期。時(shí)間步長(zhǎng)300s 實(shí)際時(shí)間步長(zhǎng)隨時(shí)間變化,由于網(wǎng)格較小,為了滿足CFL穩(wěn)定準(zhǔn)則,實(shí)際計(jì)算步長(zhǎng)在1～3s之間。水動(dòng)力參數(shù)渦流黏度0.28采用Smagorinsky方程默認(rèn)值。糙率系數(shù)32m1/3/s默認(rèn)值,介于河流糙率系數(shù)適宜范圍內(nèi)。水質(zhì)參數(shù)擴(kuò)散系數(shù)1.0同水動(dòng)力渦黏系數(shù)的1.0倍。衰減系數(shù)COD:0.2d-1TP:0.15d-1NH3-N:0.2d-1參照相關(guān)研究成果取初始值,通過(guò)模型試算求得最優(yōu)值。

根據(jù)《水域納污能力計(jì)算規(guī)程》[8]，舉水干流為污染物在橫斷面上均勻混合的中、小型河段，適宜采用河流一維模型，計(jì)算公式如下：

(1)河段的污染物濃度按式計(jì)算：

(10)

式中：C0為初始斷面的污染物濃度，mg/L；Cx為流經(jīng)x距離后的污染物濃度，mg/L；x為沿河段的縱向距離，m；u為設(shè)計(jì)流量下河段斷面的平均流速，m/s；K為污染物綜合衰減系數(shù)，s-1。

(2)相應(yīng)的水域納污能力按式計(jì)算：

M=(Cs-Cs) (Q+Qp)

(11)

式中：Q為初始斷面的入流流量，m3/s；Qp為廢污水的排放流量，m3/s；M為水域納污能力，g/s。

(3)入河排污口位于計(jì)算河段的中部時(shí)(即x=L/2)，水功能區(qū)下斷面的污染物濃度及其相應(yīng)的水域納污能力分別按以下公式計(jì)算：

(12)

M=(Cs-Cx=L)(Q+Qp)

(13)

式中：M為污染物入河速率，g/s；Cx=L為水功能區(qū)下斷面污染物濃度，mg/L。

4.2 參數(shù)確定

本研究河段為舉水干流中游段，上斷面為麻城水文站，下斷面為柳子港水文站。

(1)設(shè)計(jì)水文條件。選取2006-2016年實(shí)測(cè)水文資料，計(jì)算舉水麻城段最枯月平均流量，以作為現(xiàn)狀水平年的河道生態(tài)基流，即5.66 m3/s。根據(jù)水動(dòng)力模型分析計(jì)算，在設(shè)計(jì)水文條件下，計(jì)算河段內(nèi)平均流速均接近0.06 m/s。

(2)河流水質(zhì)目標(biāo)濃度Cs和初始斷面污染物濃度C0r。舉水水質(zhì)目標(biāo)按地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)[9]Ⅲ類計(jì)。各分區(qū)水質(zhì)目標(biāo)濃度及初始斷面污染物濃度見(jiàn)表3。

表3 舉水干流各分區(qū)特性表

Tab.3 Characteristics in all divisions

分區(qū)現(xiàn)狀水質(zhì)類別初始濃度C0r/(mg·L-1)CODNH3-NTP目標(biāo)濃度Cs/(mg·L-1)CODNH3-NTPP1Ⅱ150.150.02150.50.1P2Ⅱ150.500.10150.50.1P3Ⅱ150.500.10150.50.1P4Ⅳ150.500.10301.50.3P5Ⅳ301.500.30201.00.2P6Ⅲ301.500.30201.00.2P7Ⅲ201.000.20201.00.2P8Ⅱ201.000.20150.50.1P9Ⅱ150.500.10150.50.1P10Ⅱ150.500.10150.50.1

(3)排污口概化。計(jì)算河段匯水區(qū)內(nèi)匯集了P3～P7分區(qū)的廢污水和主要污染物。采用排污口中點(diǎn)概化的方法，得到計(jì)算河段廢污水排放流量為4.39 m3/s，COD、TP、NH3-N入河速率分別為378.27、5.38、57.53 g/s。

(4)污染物綜合衰減系數(shù)。污染物進(jìn)入水體后，水體和污染物之間同時(shí)進(jìn)行著物理、化學(xué)和生物作用的過(guò)程[10]，即污染物綜合衰減過(guò)程。污染物綜合衰減系數(shù)對(duì)計(jì)算水體的納污能力有著重要的影響。根據(jù)相關(guān)研究，我國(guó)河流COD的衰減系數(shù)為0.009～0.470 d-1，氨氮的衰減系數(shù)為0.105～0.350 d-1。根據(jù)漢江相關(guān)研究和MIKE21水動(dòng)力水質(zhì)模擬分析，綜合優(yōu)選出計(jì)算河段COD、TP和NH3-N的綜合衰減系數(shù)分別為0.20、0.15和0.20 d-1。

4.3 計(jì)算結(jié)果

(1)納污能力核算結(jié)果。根據(jù)確定的計(jì)算模型和參數(shù)，得到舉水干流各分區(qū)納污能力核算結(jié)果(見(jiàn)表4)。

(2)核算成果合理性分析。①基本資料的合理性檢驗(yàn)：采用近10年日尺度麻城水文站、柳子港水文站的實(shí)測(cè)水文資料確定納污能力核算的設(shè)計(jì)水文條件，具備合理性；水質(zhì)資料來(lái)源于湖北省環(huán)保廳地表水考核斷面監(jiān)測(cè)評(píng)估結(jié)果，具備合理性；陸域污染源資料基于當(dāng)?shù)亟?jīng)濟(jì)社會(huì)發(fā)展水平、產(chǎn)業(yè)結(jié)構(gòu)現(xiàn)狀及規(guī)劃、GDP、廢污水處理水平等資料分析得出，具備合理性。②計(jì)算條件簡(jiǎn)化和假定的合理性分析。本研究綜合考慮舉水水功能區(qū)劃和行政分區(qū)，將其分為10個(gè)分區(qū)，具備合理性。③數(shù)學(xué)模型選用、參數(shù)確定的合理性分析與檢驗(yàn)。根據(jù)《水域納污能力計(jì)算規(guī)程》，采用河流一維模型計(jì)算各分區(qū)的納污能力，水體特征污染物選COD、TP和NH3-N為指標(biāo)。污染物降解系數(shù)綜合參考相似河道的參數(shù)取值和本研究水動(dòng)力水質(zhì)模型率定結(jié)果，具備合理性。④與河流現(xiàn)狀排污情況比較分析。據(jù)調(diào)查，生產(chǎn)生活方式的不同使得各河段水質(zhì)指標(biāo)現(xiàn)狀有所差異，其差異與納污能力核算的結(jié)果相吻合。

表4 舉水主要污染物納污能力核算成果

Tab.4 Result of the water environmental capacity calculationin the Jushui River

計(jì)算分區(qū)功能區(qū)類別水質(zhì)達(dá)標(biāo)情況納污能力CODNH3-NTPP1保護(hù)區(qū)水質(zhì)達(dá)標(biāo)992.73142.5014.06P2保留區(qū)水質(zhì)達(dá)標(biāo)1378.73201.4519.57P3開(kāi)發(fā)利用區(qū)水質(zhì)達(dá)標(biāo)483.9480.266.66P4開(kāi)發(fā)利用區(qū)水質(zhì)達(dá)標(biāo)7157.30446.6586.17P5開(kāi)發(fā)利用區(qū)水質(zhì)不達(dá)標(biāo)000P6保留區(qū)水質(zhì)達(dá)標(biāo)3140.03459.1441.38P7保留區(qū)水質(zhì)達(dá)標(biāo)1282.24198.6417.67P8開(kāi)發(fā)利用區(qū)水質(zhì)達(dá)標(biāo)715.22118.7210.52P9保留區(qū)水質(zhì)達(dá)標(biāo)2330.24310.2028.81P10保留區(qū)水質(zhì)達(dá)標(biāo)1984.16308.6322.47合計(jì)19464.592266.19247.31

綜上所述，舉水納污能力核算成果具備合理性。

(3)污染物總量限排方案。通過(guò)對(duì)比分析各分區(qū)主要污染物納污能力和入河量，為達(dá)到污染物總量控制目標(biāo)，舉水流域污染物削減工作還應(yīng)進(jìn)一步加強(qiáng)。

位于麻城開(kāi)發(fā)利用區(qū)過(guò)渡區(qū)的P5分區(qū)污染物削減任務(wù)最重，與舉水干流各水功能區(qū)水質(zhì)現(xiàn)狀結(jié)果一致。該分區(qū)需以核定的納污能力作為限排總量開(kāi)展污染物削減工作，COD、NH3-N、TP的削減量分別為954.77、141.65、12.94 t/a，重點(diǎn)削減對(duì)象為工業(yè)污染源以及居民生活源。

5 結(jié) 論

本文根據(jù)舉水干流水文水質(zhì)資料，建立了MIKE21水質(zhì)模型，模擬了舉水干流麻城～柳子港河段中COD、NH3-N和TP的水質(zhì)變化過(guò)程，并按照計(jì)算規(guī)程分析計(jì)算了舉水干流各分區(qū)納污能力，結(jié)果表明：

(1)模型通過(guò)水動(dòng)力水質(zhì)耦合模擬，得到了河段內(nèi)污染物濃度變化過(guò)程圖，通過(guò)與實(shí)測(cè)水質(zhì)資料對(duì)比，模擬結(jié)果能較好反映區(qū)域內(nèi)水質(zhì)變化規(guī)律。

(2)總體上舉水干流具有一定的納污能力，COD、NH3-N和TP的納污能力分別為19 464.59、2 266.19和247.31 t/a。

(3)根據(jù)污染物限排總量分析，舉水流域亟須進(jìn)行污染物量削減工作，重點(diǎn)削減對(duì)象為工業(yè)污染源以及居民生活源。

根據(jù)舉水干流水質(zhì)現(xiàn)狀，提出以下幾點(diǎn)討論和建議：

(1)加強(qiáng)城鄉(xiāng)污水處理工程建設(shè)，盡快實(shí)施舉水河沿岸街道、鄉(xiāng)鎮(zhèn)污水集中處理全覆蓋工程。

(2)加大農(nóng)業(yè)面源治理力度，大力開(kāi)展生態(tài)濕地和生態(tài)穩(wěn)定塘等建設(shè)以提高生活廢污水污染物去除率，大力推廣生態(tài)循環(huán)農(nóng)業(yè)發(fā)展模式以提高農(nóng)田農(nóng)藥化肥利用率。

(3)流域內(nèi)麻城市、團(tuán)風(fēng)縣和新洲區(qū)存在工業(yè)聚集區(qū)，建議完善工業(yè)廢污水監(jiān)控體系，強(qiáng)化工業(yè)廢污水排放管理。

(4)建議加強(qiáng)舉水干流各水功能區(qū)交界斷面水質(zhì)監(jiān)測(cè)力度，各水功能區(qū)主管部門應(yīng)充分溝通、協(xié)調(diào)，確保交界斷面水質(zhì)達(dá)標(biāo)，為河道水質(zhì)改善提供基礎(chǔ)保障。