鄱陽湖九江工業(yè)用水區(qū)納污能力研究

楊芳 李建 裴中平 周琴

摘要：水功能區(qū)水域納污能力是污染物入河湖總量控制的基礎(chǔ)，當前水域納污能力計算較多使用公式的解析解法，計算結(jié)果較為粗略。選取鄱陽湖九江工業(yè)用水區(qū)作為研究對象，采用解析解法和數(shù)值解法相結(jié)合的方法對鄱陽湖九江工業(yè)用水區(qū)納污能力進行聯(lián)合求解，利用解析解法計算水域納污能力初值，再以此作為模型啟動條件，通過數(shù)值模型對解析解法獲得的初始值進行試算檢驗和優(yōu)化。計算結(jié)果表明，鄱陽湖九江工業(yè)用水區(qū)的最優(yōu)水域納污能力COD為9 776 t/a，氨氮為142 t/a。研究統(tǒng)籌考慮了河道形態(tài)、岸邊流速變化、豐水期長江江水頂托和倒灌等因素，彌補了解析解法的局限性，結(jié)果更具合理性，可為鄱陽湖九江工業(yè)用水區(qū)水質(zhì)管理和限制排污總量控制提供參考。

關(guān)鍵詞：水域納污能力;解析解法;數(shù)值解法;鄱陽湖;工業(yè)用水區(qū)

中圖分類號：X52文獻標志碼：A

開放科學（資源服務(wù)）標識碼（OSID）：李建

Assimilative capacity of Jiujiang industrial water function zone beside Poyang Lake using analytical method and numerical method

YANG Fang，LI Jian，PEI Zhongping，ZHOU Qin

（Changjiang Water Resources Protection Institute，Wuhan 430051，China）

Abstract：The assimilative capacity of a water function zone is the basis for gross control of pollutants discharge.At present，the capacity is commonly calculated by analytical method，but the calculation results are not accurate.In this research，we took Jiujiang industrial water function zone beside Poyang Lake as the object of research，and calculated its assimilative capacity.The analytical method was used to calculate the initial value，and the numerical method was used to optimize the initial value.The calculation results showed that the optimal assimilative capacity of Jiujiang industrial water function zone beside Poyang Lake for COD is 9 776 t/a，and that for ammonia nitrogen （NH3-N） is 142 t/a.This study takes into account such factors as the river morphology，the change of inshore velocity，and the effects of characteristic hydrological conditions.It can make up for the limitations of the analytical method and generate more reasonable results.It can provide a reference for water pollutants gross control in Jiujiang industrial water function zone beside Poyang Lake.

Key words：assimilative capacity of water bodies;analytical method;numerical method;Poyang Lake;industrial water function zone

水域納污能力計算是江河湖泊限制排污總量控制的基礎(chǔ)。《中國水利百科全書》（2004年）環(huán)境水利分冊將水域納污能力定義為水域能持續(xù)發(fā)揮給定功能而接納的最大污染物負荷量，即根據(jù)規(guī)定的水文概率或經(jīng)調(diào)節(jié)的河道流量，確定在自凈作用下，保證達到功能區(qū)水質(zhì)要求所允許的污染物排放量。《水域納污能力計算規(guī)程》（GB/T 25173-2010）明確水域納污能力是指，針對江河、湖泊、水庫、運河、渠道等已劃定水功能區(qū)的地表水域，在設(shè)計水文條件下，某種污染物滿足水功能區(qū)水質(zhì)目標要求所能容納的該污染物的最大數(shù)量。

國外常用“環(huán)境容量”、“最大容許納污量”和“水體容許排污水平”等概念來描述水域納污能力，并將其融合在總量控制研究中[1]。20世紀60年代，日本最早開始水體允許污染物負荷量估算研究，提出“環(huán)境容量”概念[2]。隨后歐美國家相繼開展環(huán)境容量研究，通常將水環(huán)境容量計算和污染物總量分配在同一過程中實現(xiàn)[3]，早期較多采用公式法、線性規(guī)劃等確定性方法進行研究[4-5]，未考慮河流水文、水動力及水質(zhì)條件的動態(tài)變化，之后隨機理論被引入水環(huán)境容量研究，先后提出了概率稀釋模型、隨機水質(zhì)優(yōu)化模型等不確定方法[6-8]，當前研究較多關(guān)注于河流流量、排污流量、排污濃度、河流水質(zhì)背景濃度等因子不確定性對污染負荷分配的影響。

國內(nèi)對于流域水域納污能力的研究經(jīng)歷了由濃度控制到目標總量控制再到容量總量控制的發(fā)展歷程，例如，太湖流域已將容量總量控制運用于重點控制區(qū)的污染防治。目前針對水域納污能力的研究方法較多，國內(nèi)學者改進了概率稀釋模型法[9]，提出了模型試錯法[10-12]和基于盲數(shù)理論的未確知數(shù)學法[13]等計算方法，但在管理實踐中仍以公式法[14]、線性規(guī)劃法[15]等確定性方法為主，其主要原因是國內(nèi)水環(huán)境管理體制、水污染現(xiàn)狀與國外存在較大差異，且多數(shù)地區(qū)缺乏長系列水質(zhì)監(jiān)測數(shù)據(jù)[16-17]，國外先進方法在國內(nèi)的適用條件尚不夠成熟。

3計算步驟與結(jié)果

綜合利用兩種方法進行水域納污能力的計算，通過解析解法求得水功能區(qū)納污能力的初始值，再以此作為數(shù)值模型起算條件，利用數(shù)值解法對初始值進行試算檢驗和優(yōu)化。

3.1設(shè)計水文條件與計算參數(shù)

結(jié)合鄱陽湖“高水是湖、低水似河，洪水一片、枯水一線”的鮮明特點，采用湖口水文站90%最枯月均流量作為設(shè)計水文條件。由于鄱陽湖在汛期存在江水倒灌情況，在統(tǒng)計最枯月均流量時排除了江水倒灌情況下的最枯月流量。對1950-2015年湖口水文站流量資料進行統(tǒng)計，用適線法獲得湖口水文站90%最枯月均流量為632 m.3/s，對應(yīng)的多年平均水位為6.76 m（吳淞高程）。

蛤蟆石斷面位于鄱陽湖九江工業(yè)用水區(qū)內(nèi)，距離功能區(qū)上游約500 m，功能區(qū)上游斷面至蛤蟆石間沒有排污口，因此選取蛤蟆石斷面水質(zhì)作為功能區(qū)背景濃度控制斷面。根據(jù)近年（2014-2017年）蛤蟆石斷面的水質(zhì)監(jiān)測結(jié)果，枯水期（11月至次年2月）COD平均濃度10.17 mg/L，氨氮0.46 mg/L，以此作為數(shù)值模型模擬的背景濃度。

盡管鄱陽湖九江工業(yè)用水區(qū)水質(zhì)管理目標為Ⅳ類，但其下游水功能區(qū)水域——鄱陽湖湖區(qū)保留區(qū)的水功能區(qū)水質(zhì)目標為Ⅱ類，因此鄱陽湖九江工業(yè)用水區(qū)下斷面必須以Ⅱ類水質(zhì)目標嚴格控制，即COD取15 mg/L，氨氮取0.5 mg/L。

3.2基于解析解法的納污能力初值計算

根據(jù)水域納污能力解析解計算公式（8），得到在設(shè)計水文條件下鄱陽湖九江工業(yè)用水區(qū)水域納污能力的初始值，COD為8217.6 ?t/a，氨氮為85.25 t/a（見表1）。

3.3基于數(shù)值解法的納污能力終值試算

3.3.1數(shù)值模型的建立與驗證

采用MIKE21模型進行數(shù)值計算。鄱陽湖九江工業(yè)用水區(qū)下游為湖口水文站，上游為星子水位站，模型模擬范圍選定為星子-湖口段，長40 km。

采用三角形結(jié)合四邊形網(wǎng)格的方式對鄱陽湖星子-湖口段40 km計算水域（1∶10000水下地形，鄱陽湖水文局提供）[HJ2.1mm]進行剖分。其中，鄱陽湖九江工業(yè)用水區(qū)河段順直，采用四邊形網(wǎng)格，其他非順直河段采用三角形網(wǎng)格以適應(yīng)復(fù)雜的湖泊岸線邊界。四邊形網(wǎng)格長50～60 m，寬20～30 m，三角形最大面積10 000 m.2，采用鄰近點法對網(wǎng)格節(jié)點高程進行插值，形成水動力水質(zhì)模型的計算網(wǎng)格（圖2）。

水動力模型和對流擴散模型參數(shù)采用經(jīng)驗公式、以往研究成果和模型試算等多種方式綜合確定，各參數(shù)取值見表2所示。

采用2015年屏峰站水位和蛤蟆石斷面水質(zhì)實測數(shù)據(jù)對模型進行驗證，數(shù)值模型的驗證結(jié)果顯示模擬誤差在可接受范圍內(nèi)（圖3）。

3.3.2污水排放量固定時允許排放濃度的試算[HJ1.9mm]

為了便于計算，通常需要將水功能區(qū)排污口概化至河段中間位置。在水域納污能力初始值M0確定的情況下，污水排放量Q污與污染物濃度C污有

多重反比例組合關(guān)系。在利用模型數(shù)值解法進行污染物擴散情況計算的時候，首先需要確定一個起算的污水流量值，由此對應(yīng)一個污染物的起算濃度。在概化排污口一系列流量Q污與濃度C污組合中，篩選C污值的條件是判斷概化排污口排污后污染物到達水功能區(qū)下斷面時，其濃度是否剛好達到水質(zhì)控制目標的上限。

根據(jù)污染源現(xiàn)狀調(diào)查成果和近期規(guī)劃，研究區(qū)域近期入湖污水總量將達到22.77萬t/d，即污水排放流量可達2.64 m.3/s。在預(yù)留部分發(fā)展空間基礎(chǔ)上，確定鄱陽湖九江工業(yè)用水區(qū)納污能力試算的污水排放量初值為3 m.3/s，由此可推出，概化排污口COD的對應(yīng)濃度約為87 mg/L，氨氮對應(yīng)濃度約為0.9 mg/L，以此確定不同濃度的試算工況見表3。

對多組污染物輸入濃度進行試算，取各代表斷面污染物最大濃度統(tǒng)計不同污染物輸入濃度工況下COD濃度、氨氮濃度沿程變化（表4）。

從表4中可以看出，當概化排污口污水排放流量為3 m.3/s，COD排放濃度為95 mg/L，氨氮排放濃度為1.5 mg/L時，鄱陽湖九江工業(yè)用水區(qū)下斷面、距離概化排污口下游2 500 m處COD和氨氮濃度均可恢復(fù)至地表水Ⅱ類控制目標（分別為15 mg/L和0.5 mg/L）。在上述污水量和污染物輸入條件下，計算得出概化排污口的污染物入湖負荷量即為鄱陽湖九江工業(yè)用水區(qū)的水域納污能力，其中COD為8 988 t/a，氨氮為142 t/a。

3.3.3污染物排放濃度固定時允許排放流量的試算

為嚴格控制入河污染物濃度，鄱陽湖九江工業(yè)用水區(qū)要求現(xiàn)有排污口污水排放標準執(zhí)行《城鎮(zhèn)污水處理廠污染物排放標準》一級A標準，即COD排放濃度為50 mg/L，氨氮排放濃度為5 mg/L。上述試算計算結(jié)果中，氨氮已滿足污水排放標準要求，在該排污濃度下計算得出的氨氮納污能力數(shù)值已經(jīng)是一個嚴格限值。而COD為 95 mg/L的排放濃度不滿足污水排放標準，因此需繼續(xù)優(yōu)化COD的排放濃度和排放量。

表4計算結(jié)果顯示，在污水排放流量一定的情況下，污染物排放濃度越小，對排污口下游水質(zhì)的影響越小，當COD排放濃度為95 mg/L時，鄱陽湖九江工業(yè)用水區(qū)下斷面（概化排污口下游2 500 m）處

COD濃度已衰減至15 mg/L，由此可推斷COD排放濃度執(zhí)行《城鎮(zhèn)污水處理廠污染物排放標準》一級

A標（COD≤50 mg/L）時，概化排污口下游2 500 m處COD濃度低于15 mg/L，優(yōu)于Ⅱ類水質(zhì)管理目標，此時該水功能區(qū)的承載能力沒有得到完全利用。基于此，以COD排放濃度50 mg/L為模型輸入條件，以鄱陽湖九江開發(fā)利用區(qū)下斷面COD≤15 mg/L為控制條件，試算可允許最大污水排放量。

通過對多組流量輸入條件進行試算，取各代表斷面污染物濃度最大值統(tǒng)計得到不同工況下COD濃度沿程變化表，見表5。

試算結(jié)果表明，當概化排污口COD濃度為50 mg/L，污水流量為6.2 m.3/s時，鄱陽湖九江工業(yè)用水區(qū)下斷面COD可滿足Ⅱ類水質(zhì)控制目標（15 mg/L）。在上述污染物排放濃度和污水排放量條件下，概化排污口COD入湖負荷量即為鄱陽湖九江工業(yè)用水區(qū)COD的納污能力，為9 776 t/a。

綜合上述試算結(jié)果，得到鄱陽湖九江工業(yè)用水區(qū)的最優(yōu)水域納污能力COD為 9 776 t/a，氨氮為142 t/a。污染物按該納污能力負荷量排放時COD和氨氮影響范圍，見圖4。

4討論

4.1解析解法的局限性分析

解析解法是目前普遍采用的一種開發(fā)利用區(qū)水域納污能力計算方法，也是《水域納污能力計算規(guī)程》（GB/T 25173-2010）推薦使用的一種方法。但是，解析解法也存在一定的局限性[26]。

（1）河流二維模型橫斷面概化的局限性。進行解析解計算時必須假定污染物連續(xù)恒定排放，計算河段橫斷面為矩形，對于橫斷面不能概化為矩形的河段，原則上只能通過數(shù)值模型求解計算。鄱陽湖湖口水道地形復(fù)雜，復(fù)式斷面特征明顯，河道岸灘地形變化造成的影響，解析解法無法精確反應(yīng)。利用二維模型進行解析解求解的另一重要特點條件是假定岸邊計算水域斷面流速為恒定均勻的，這種概化忽略了斷面橫向的流速變化，因此通過解析解獲得的水域納污能力值是較為粗略的。

（2）計算水域污染帶寬度取值的局限性。對于[CM（22]河寬比較大的大型河流，污染物從岸邊排放后不可能達到全斷面混合，即容許污染物降解稀釋的混合區(qū)一般在岸邊一定范圍內(nèi)，此時需要計算岸邊水域的納污能力。在不使用數(shù)值模型模擬的情況下，岸邊水域計算范圍的取值存在較大不確定性，同樣對水域納污能力計算結(jié)果的準確性產(chǎn)生一定影響。

解析解法獲得的納污能力比較適合作為數(shù)值解模型試算的起算條件，是納污能力精確計算的重要參考，有助于快速確定水功能區(qū)納污能力的最優(yōu)值。

4.2豐水期壅水條件影響分析

鄱陽湖是長江中下游最大的通江湖泊，與長江之間存在著復(fù)雜的水文和水動力交互作用，而江水倒灌是長江頂托過程的極端現(xiàn)象。鄱陽湖壅水后，湖口和星子之間的水位差變小，流速減小，對污染物遷移擴散可能產(chǎn)生不利影響，因此有必要進一步討論特征豐水期水文條件下的排污影響，即保持概化排污口排污條件（最佳排放流量和濃度）不變，分析特征水動力條件下的污染物擴散影響。

一般情況下，當湖口與星子之間的水位差為0時，鄱陽湖九江工業(yè)用水區(qū)水流條件對污染物遷移擴散較為不利。當湖口與星子之間水位差為0時有兩種情況發(fā)生，一種是長江高水位對鄱陽湖水面產(chǎn)生頂托，但是沒有發(fā)生江水倒灌現(xiàn)象;另一種是在江水頂托基礎(chǔ)上進一步發(fā)生了江水倒灌現(xiàn)象。篩選1950-2015年期間湖口與星子之間水位差為0的日期，發(fā)現(xiàn)1971年7月30日和1963年8月25日最具代表性，此時鄱陽湖水位處于相對較低的狀態(tài)，研究區(qū)域的水動力條件也相對不利，故以此作為豐水期的工況條件（表6）。

計算豐水期兩組工況，統(tǒng)計分析排污口下游沿程污染物最大濃度變化，結(jié)果見圖5和表7。模擬計算表明，豐水期在長江江水頂托、倒灌等特征水文條件下，按照COD 9 776 t/a、氨氮142 t/a污染物入湖量計算，鄱陽湖開發(fā)利用區(qū)下斷面COD和氨氮濃度均優(yōu)于水質(zhì)控制目標。因此，本研究計算的COD和氨氮納污能力在豐水期特征水文條件下也是符合水質(zhì)控制條件的。

5結(jié)論

國內(nèi)水功能區(qū)管理實踐中，開發(fā)利用區(qū)納污能力計算大多采用公式的解析解法，但在模型概化、污染帶寬度取值等方面存在較多局限性。本文利用公式的解析解法與數(shù)值解法相結(jié)合的方法對鄱陽湖九江開發(fā)利用區(qū)納污能力進行了聯(lián)合求解。通過研究主要得出以下結(jié)論。

（1）平面二維水動力水質(zhì)方程的求解方法分為解析解法和數(shù)值解法，其中解析解法計算結(jié)果較為粗略，而數(shù)值解法計算結(jié)果更為準確。利用解析解法計算結(jié)果作為數(shù)值解法試算的起算條件，有助于快速確定水功能區(qū)納污能力的最優(yōu)值。

（2）利用數(shù)值解法進行計算時，在污染物排放總負荷量一定情況下，污水排放量和污染物排放濃度存在反比例組合關(guān)系，需以水功能區(qū)下斷面水質(zhì)管理需求作為判斷條件，遴選污水排放量和污染物排放濃度的最佳組合，同時需綜合考慮河道形態(tài)、岸邊水深和流速變化、特征水文條件影響，以及污染物排放標準要求等條件。

（3）經(jīng)計算，鄱陽湖九江開發(fā)利用區(qū)納污能力最優(yōu)值為COD為9 776 t/a，氨氮為142 t/a。

參考文獻（References）：

[1]羅小勇主編.基于水功能區(qū)的納污能力計算理論方法及應(yīng)用[M].北京：中國水利水電出版社，2016.（LUO X Y.Theoretical method and application of allowable permitted Assimilative Capacity based on water function zone[M].Beijing：China Water & Power Press，2016.（in Chinese））

[2]夏青.水環(huán)境容量[J].環(huán)境保護科學，1981（4）：21-29.（XIA Q.Water environment capacity[J].Environmental Protection Science，1981（4）：21-29.（in Chinese）） DOI：10.16803/j.cnki.issn.1004-6216.1981.04.005.

[3]董飛，劉曉波，彭文啟，等.地表水水環(huán)境容量計算方法回顧與展望[J].水科學進展，2014，25（3）：451-463.（DONG F，LIU X B，PENG W Q，et al.[JP+4]Calculation methods of water environmental capacity of surface waters：Review and prospect[J].Advances in Water Science，2014，25（3）：451-463.（in Chinese）） DOI：10.14042/j.cnki.32.1309.2014.03.020.

[4]LIEBMAN J C，LYNN W R.The optimal allocation of stream dissolved oxygen[J].Water Resource Research，1966，2（3）：581-591.DOI：10.1029/WR002i003p00581.

[5]REVERLLE C S，LOUCKS D P，LYNN W R.Linear programming applied to water quality management[J].Water Resource Research，1968，4（1）：1-9.DOI：10.1029/WR004i001p00001.

[6]FUJIWARA O，GNANENDRAN S K，OHGAKI S.River quality management under stochastic streamflow[J].Journal of Environmental Engineering，1986，112（2）：185-198.DOI：10.1061/（ASCE）0733-9372（1986）112：2（185）.

[7]USPA.Technical guidance manual for performing wasteload allocation，book：VII：Permit averaging periods[R].Washington DC：United States Environmental Protection Agency，Office of Water，1984：1-4.

[8]ELLIS J H，MCBEAN E A.Optimization modeling of water quality in an uncertain environment[J].Water Resource Research，1987，23（12）：2227-2238.DOI：10.1029/WR021i007p00934.

[9]夏青，孫艷，賀珍，等.水污染物總量控制實用計算方法概要[J].環(huán)境科學研究，1989，2（3）：1-73.（XIA Q，SUN Y，HE Z.Summary of practical calculation methods for total amount control of water pollutants[J].Research of Environmental Sciences，1989，2（3）：1-73.（in Chinese）） DOI：10.13198/j.res.1989.03.3.xiaq.001.

[10][ZK（#]韓能進.河流一維水質(zhì)模型在水環(huán)境容量計算方面的應(yīng)用[J].環(huán)境科學與技術(shù)，1995（4）：43-45.（HAN N J.Application of one dimensional river water quality model to water environmental capacity calculation[J].Environmental Science and Technology，1995（4）：43-45.（in Chinese）） DOI：10.19672/j.cnki.1003-6504.1995.04.013

[11]曾思育，徐一劍，張?zhí)熘?環(huán)狀河網(wǎng)水質(zhì)模型在水污染控制規(guī)劃中的應(yīng)用[J].水科學進展，2004，15（2）：193-196.（ZENG S Y，XU Y J，ZHAGN T Z.Application of unsteady water quality model for looping river network to water pollution control planning[J].Advances in Water Science，2004，15（2）：193-196.（in Chinese）） DOI：10.14042/j.cnki.32.1309.2004.02.012.

[12]張紅舉，甘升偉，袁洪州，等.環(huán)太湖河流入湖水質(zhì)控制濃度分析[J].水資源保護，2012，28（6）：8-11.（ZHANG H J，GAN S W，YANG H Z，et al.Analysis of controlled concentration of water quality of rivers around Taihu Lake[J].Water Resources Protection，2012，28（6）：8-11.（in Chinese）） DOI：10.3969/j.issn.1004-6933.2012.05.002.

[13]李如忠，汪家權(quán)，王超，等.不確定性信息下的河流納污能力計算初探[J].水科學進展，2003，14（4）：359-363.（LI R Z，WANG J Q，WANG C，et al.Calculation of river water environmental capacity under unascertained information[J].Advances in Water Science，2003，14（4）：359-363.（in Chinese）） DOI：10.14042/j.cnki.32.1309.2003.04.013.

[14]周娜，賈仰文，胡鵬，等.松花江流域冰封期水功能區(qū)限制納污能力控制研究[J].水利學報，2014，45（5）：557-565.（ZHOU N，JIA Y W，HU P.Research on the regulation of water function reaches pollution load in freeze up period of Songhuajiang River Basin[J].Journal of Hydraulic Engineering，2014，45（5）：557-565.（in Chinese）） DOI：10.13243/j.cnki.slxb.2014.05.007.

[15]張守平，魏傳江，王浩，等.流域/區(qū)域水量水質(zhì)聯(lián)合配置研究I：理論方法[J].水利學報，2014，45（7）：757-766.（ZHANG S P，WER C J，WANG H，et al.Basin/region water quality allocation I.Theory and method[J].Journal of Hydraulic Engineering，2014，45（7）：757-766.（in Chinese）） DOI：10.13243/j.cnki.slxb.2014.07.001.

[16]張曉，羅軍剛，解建倉.考慮取水口和支流的河流納污能力計算模型研究與應(yīng)用[J].水利學報，2017，48（3）：317-324.（ZHANG X，LUO J G，XIE J C.Study on the calculation model of river water environment capacity considering water intake and tributary[J].Journal of Hydraulic Engineering，2017，48（3）：317-324.（in Chinese）） DOI：10.13243/j.cnki.slxb.20160270.

[17]王新娟，孫穎，邵景力，等.北京西郊地區(qū)地下水納污能力分析[J].南水北調(diào)與水利科技，2014，12（6）：31-34.（WANG X J，SUN Y，SHAO J L，et al.Pollutant carrying capacity analysis of groundwater in the Western Suburbs of Beijing[J].South-to-North Water Transfers and Water Science & Technology，2014，12（6）：31-34.（in Chinese）） DOI：10.13476/j.cnki.nsbdqk.2014.06.007.

[18]芮孝芳.水文學原理[M].北京：中國水利水電出版，2004.（RUI X F.Physical Hydrology[M].Beijing：China Water & Power Press，2004.（in Chinese））

[19]裴中平，楊芳，辛小康.入河排污口設(shè)置論證技術(shù)與實例[M].武漢：長江出版社，2017.（PEI Z P，YANG F，XIN X K.Technology and practice For demonstration on outfall setting up[M].Wuhan：Changjiang Press，2017.（in Chinese））

[20]逄勇，陸桂華.水環(huán)境容量計算理論及應(yīng)用[M].北京：科學出版社，2016.（PANG Y，LU G H.Theory and application of water environmental capacity calculation[M].Beijing：Science Press，2016.（in Chinese））

[21]蘇麗娟.雙邊空間分數(shù)階對流擴散方程的幾種數(shù)值解法[D].濟南：山東大學，2010.（SU L J.Some numerical methods for fractional diffusion equations[D].Jinan：Shangdong University，2010.（in Chinese））

[22]謝蘭英，奚紅霞，謝新平，等.不同邊界條件下球形顆粒吸附劑內(nèi)擴散方程的分析解和數(shù)值解的比較[J].離子交換與吸附，2000，16（4）：289-295.（XIE L Y，XI H X，XIE X P，et al.A comparison between analytical solutions and numerical solution of diffusion equation in spherical adsorbent particle under different boundary conditions[J].Ion Exchange and Adsorption，2000，16（4）：289-295.（in Chinese）） DOI：10.16026/j.cnki.iea.2000.04.001.

[23]GB/T 25173-2010.水域納污能力計算規(guī)程[S].2010.（GB/T 25173-2010.Code of practice for computation on allowable permitted assimilative capacity of water bodies[S].2010.（in Chinese））

[24]ZHAO D H，SHEN H W，TABIOS G Q，et al.Finite-volume two-dimensional unsteady-flow model for river basins[J].Hydr Engrg，ASCE.1994，120（7）：863-883.DOI：10.1061/（ASCE）0733-9429（1994）120：7（863）

[25]SLEIGH P A，GASKELL P H，BERZINS M，et al An unstructured finite-volume algorithm for predicting flow in rivers and estuaries[J].Computers & Fluids，1998，27（4）：479-508.DOI：10.1016/S0045-7930（97）00071-6.

[26]馬莉，桂和榮，曹彭強.河流污染二維水質(zhì)模型研究及RMA4模型概述[J].安徽大學學報（自然科學版），2011，35（1）：102-108.（MA L，GUI H R，CAO P Q.Study on integrateion of 2D water quality models and review of RM A4 model[J].Journal of Anhui University （Natural Science Edition），2011，35（1）：102-108.（in Chinese）） DOI：10.3969/j.issn.1000-2162.2011.01.020.