基于不同環(huán)境下南渡江河口污染物分布研究

方明儀，施浩然*，海來以波，李海生，明 睿

(1.西華大學能源與動力工程學院，四川 成都 610039；2.西華大學流體及動力機械教育部重點實驗室，四川 成都 610039)

21世紀是“海洋的世紀”。伴隨著海洋經(jīng)濟的開發(fā)，沿海環(huán)境已經(jīng)成為社會各界關注的熱點。隨著經(jīng)濟社會的快速發(fā)展，海洋開發(fā)規(guī)模與強度持續(xù)加大，近岸海域污染日益嚴重，海域環(huán)境質量日益惡化[1]。河口水環(huán)境問題突出，沿海水域的污染壓力主要來自上游淡水攜帶污染物匯入河口。近年來，計算機技術的發(fā)展，提高了對數(shù)據(jù)的處理效率，水質模型也用得越來越頻繁。國內對河口及突發(fā)事件水質問題的研究較多，趙棣華等[2]采用潮流水質模型模擬了長江江蘇省感潮河段COD的濃度分布，得出模擬得到的COD濃度范圍與圖片遙感分析結果基本一致；徐帥[3]、李大鳴[4]、楊晨等[5]基于MIKE水動力-水質模型對河流及水庫污染物在水體中對流擴散進行了研究，得出在不同工況下各類污染物的影響及分布范圍；陳冬云等[6]采用潮流水質模型研究了徑流量對錢塘江河口突發(fā)污染物遷移擴散的作用，結果表明，加大上游徑流量可明顯加快污染團的遷移、擴散和稀釋速度。

目前來說，中國對入海中小河流的基礎研究還很弱，與長江、黃河、珠江等大江大河相比，中小河流入海物質的通量變化及其對流域-河口環(huán)境效應，以及自然與人類活動雙重因素驅動下陸源物質從流域到河口的源-潮匯過程的研究不足[7]。對熱帶河口污染物擴散變化研究相對較少，尤其是像南渡江河口，受東部臺風影響較為明顯，是多臺風登陸的重風害地區(qū)，以及常年水體溫度較高，對水體中各種污染物之間的理化生反應劇烈，加上南渡江河口隱藏排放較多，河口兩岸面源排放大，導致污染源識別過程存在著很多不確定性因素。目前對于觀測數(shù)據(jù)、模型參數(shù)的不確定性方面有一定的研究成果，對污染源信息(如排放方式、排放位置)、水動力條件(如潮汐水域)、污染物性質的不確定性(如有機污染物)方面的研究成果缺乏[8]。

根據(jù)2018年海南省海洋環(huán)境狀況公報可知，在2018年有9個臺風進入南?；蛟谀虾Ｉ桑渲?，1804號臺風(熱帶風暴級)“艾云尼”和1809號臺風(熱帶風暴級)“山神”登陸海南島，在海南島沿岸引發(fā)不同程度的風暴增水，根據(jù)潮位站資料顯示，?？谛阌Ⅱ灣闭咀畲笤鏊?.38 m，瓊海博鰲驗潮站最大增水0.64 m。根據(jù)《海南省水資源公報》可知，海口市用水量加大，主要是?？诔ｑv人口逐年遞增，工農(nóng)業(yè)用水量加大，上游需水量逐年遞增。隨著?？谑谐鞘谢M程的不斷加快，河道面積縮減、城市排水不暢、水質惡化等水環(huán)境問題日趨嚴重[9]。

為此，根據(jù)南渡江河口水環(huán)境情況，對南渡江河口水質變化規(guī)律進行研究?；诔毕蜕嫌蝸砹飨嗷プ饔孟拢_展在不同環(huán)境下污染物在河道中變化規(guī)律。實現(xiàn)在熱帶復雜河口基礎上，探究極端天氣和惡劣條件下污染物變化情況，提出水生態(tài)保護與綜合治理建議和方案，為海南省“十四五”規(guī)劃水環(huán)境治理提供依據(jù)。

1 研究區(qū)域概況

選取黑山村斷面至北干流(南渡江下游)入?？跒檠芯糠秶?，河長16 km，河道平均坡降0.072%，研究范圍見圖1。入海口河流水系結構復雜，河道通過多處分汊河道與外海聯(lián)通，并形成多處城市水域。在入?？谔幏帚鉃?條支流，主干為北干流，受東北季風而迫使分汊向西，依次成橫溝河、海甸溪。北干流主要以排水排沙為主，排沙約占6/7，排水約占9/10；橫溝河排沙1/7，排水1/10，河長5.3 km；海甸溪河長2.8 km，是進水進沙的潮汐河道，漲潮流向東，落潮流向西[10]。河口段受潮汐作用影響，潮流界可上行至鐵橋位置，潮區(qū)界可達更遠些。

南渡江入?？诔Ｄ昶骄w溫度在25℃，有著獨特的臺風季節(jié)特征，南渡江流域降雨量分配非常不均，每年5—10月降雨量就達到全年80%，河道水量主要靠降雨補給。

2 研究方法

2.1 MIEK21模型基本原理

隨著科學技術的進步及相關理論的不斷完善，水動力模擬運用越來越廣泛，在實踐中發(fā)揮著重要作用[11]。MIKE21作為眾多模擬軟件之一，可用于模擬湖泊、河流、泥沙、海灣、海洋的波浪及水流。其水動力模塊(HD)遵循N-S方程[12]，并服從Boussinesq假設，即流體低速流動時僅考慮溫度對流體密度的影響，同時服從靜水壓力。MIKE21水動力模型的控制方程包括連續(xù)方程和動量方程，對流擴散的控制方程是對流擴散方程。

圖1 南渡江河口研究范圍

2.1.1控制方程

對于水平尺度遠大于垂直尺度的情況，由于水深、流速等水力參數(shù)沿垂直方向的變化比沿水平方向的變化要小，因此將三維流動的控制方程沿水深積分，并取水深平均，可得到沿水深平均的二維淺水流動質量和動量守恒控制方程組。其連續(xù)性方程、X和Y方向動量方程，可分別表示為：

(1)

(2)

(3)

式中h——水深，h=d+ζ；ζ、d——水位、水深，m；p、q——x、y方向上的流通通量，即單寬流量，m2/s；C——謝才系數(shù)，m1/2/s；g——重力加速度，m/s2；Ω——科氏力系數(shù)(Ω=0.729×10-4)；ρ——水的密度，kg/m3；V、Vx、Vy——風速及在x、y方向上的分量，m/s；f——風阻力系數(shù)；τxx、τxy、τyy——剪切應力分量，kg/m2；pa——當?shù)卮髿鈮簭?，Pa。

2.1.2對流擴散方程

(4)

式中c——水質濃度，mg/L；Dx、Dy——x、y方向的擴散系數(shù)，m2/s；u、v——x、y方向上的流速分量。

3 模型構建

3.1 網(wǎng)格剖分

網(wǎng)格劃分尤為重要，在MIKE21網(wǎng)格劃分工具中，網(wǎng)格劃分的好壞直接影響模型穩(wěn)定及誤差大小[13]。南渡江河口研究區(qū)，地理坐標范圍：110°16′～110°24′E，19°57′～20°05′N。由于河道曲率半徑變化不同，岸線曲折不規(guī)則，所以采用非結構性網(wǎng)格對研究區(qū)域進行三角形剖分。

網(wǎng)格共有11 157個，節(jié)點數(shù)6 220個，最小單元格面積223 m2，最大單元格面積5 550 m2。三角形網(wǎng)格最小角度25.4°。最終網(wǎng)格剖分見圖2。

圖2 三角形網(wǎng)格剖分

3.2 地形插值

利用項目區(qū)河口段1∶100 00實測水下地形圖，提取相應河道高程點，基于非結構性網(wǎng)格基礎上，導入河道高程文件，采用自然鄰近插值法進行，得到河底地形(圖3)。

圖3 河床高程分布(m)

3.3 邊界條件及源匯設置

a)初始條件。采用模型熱啟動方式，即首先對模型進行模擬，模型達到穩(wěn)定后提取模擬結果，再將模擬結果導入初始文件中作為模型的初始條件。

b)邊界。模型上邊界在黑山村斷面處，下邊界分別在南渡江、橫溝河、海甸溪入?？?。水文邊界條件的設置：上游采用龍?zhí)了恼緦崪y2016年流量數(shù)據(jù)，下游采用2016年潮位作為下游邊界；水質邊界設置：上游來流污染物濃度采用逐月平均實測生成的時間序列作為輸入條件。

c)點源。在河流分汊處，排污口通過明渠形式排入南渡江，經(jīng)過資料收集查找并分析得出在該排污口的排放量為30 m3/s，排放的污染物有總磷、氨氮、化學需氧量、高錳酸鹽，它們年均每天排放含量分別為0.08、0.52、14.04、5.03 mg/L，將此排污口概化為點源。

模型上邊界為南渡江上游，下邊界為南渡江下游、橫溝河下游、海甸溪下游，以及排污口位置見圖4。

圖4 邊界及排污口位置示意(m)

3.4 重要參數(shù)設置

3.4.1糙率

通常情況下，河道斷面平均水深較小，糙率變化較大，尤其水深較小糙率逐漸變大，當水深到達一定深度，此時糙率值最小，之后水深增加糙率幾乎不變[14]。在工程項目中，通常有3種方法確定各計算河段的糙率，包括利用河道實測水文資料推算糙率、查表法及糙率公式[15]。

對南渡江研究范圍河道現(xiàn)狀及兩岸植被分布情況進行現(xiàn)場查勘及洪水調查結果，首先根據(jù)水深與糙率的反比關系，通過線性差值，得到初始河道糙率場，再通過反復率定和驗證，確定研究范圍糙率分布(圖5)。

圖5 糙率場分布(m)

3.4.2擴散系數(shù)

影響擴散系數(shù)的主要因素是：河寬和水深。在相同水深下，河道越寬，擴散系數(shù)越大；在相同河寬下，水深越深，擴散系數(shù)越大。MIKE 21中的二維對流擴散方程假定垂向充分混合，不考慮垂向作用僅考慮水平方向的對流擴散。經(jīng)過多次率定及驗證之后，污染物擴散系數(shù)取0.15 m/s2。

3.4.3降解系數(shù)

此次降解系數(shù)采用綜合降解系數(shù)，南渡江河口段污染物降解基本符合一級降解動力反應方程[16]。計算公式為：

Ct=Coe-kt

(5)

(6)

式中t——反應時間，d；k——綜合降解系數(shù)，1/d；Ct——污染物濃度，mg/L；Co——污染物初始濃度。

根據(jù)該公式初估各種污染物的降解系數(shù)，再經(jīng)過反復驗證后，最終確定南渡江河口區(qū)KMnO4降解系數(shù)為1.5E-006/s，COD降解系數(shù)為6.5E-007/s，NH3-N降解系數(shù)為3.8E-007/s，TP降解系數(shù)為1.5E-006/s，TN降解系數(shù)為3E-007/s。

3.5 MIKE 21 HD 模擬驗證

模型水動力計算時間選取在夏季豐水時期，時間從2016年6月28日19:00至2016年6月29日23:00。通過入?？诜帚恻cL1實測值，對模型的流速、流向和水位進行驗證。流速、水位及流向(驗證)見圖6。

a)流速

根據(jù)水動力率定情況，進一步通過相關系數(shù)(R)和納什效率系數(shù)(NSE)對模擬結果進行了統(tǒng)計分析，統(tǒng)計結果見表1，流速、流向及水位的相關系數(shù)均大于0.76，納什效率系數(shù)也在0.58以上。因此，本次構建的南渡江河口區(qū)MIKE21水動力模擬結果良好，模型參數(shù)設置合理。相關系數(shù)和納什效率系數(shù)計算公式如下：

(7)

(8)

表1 水動力變量統(tǒng)計分析

3.6 MIKE 21 AD模擬及驗證

基于水動力構建基礎上，加入對流擴散模塊對4種污染物進行模擬，模擬時段與水動力一致。通過儒房水質監(jiān)測站實測資料對污染物進行驗證，通過相關系數(shù)(R)和納什效率系數(shù)(NSE)對4種污染物模擬結果進行了統(tǒng)計分析，統(tǒng)計結果見表2，高錳酸鹽、氨氮、總磷及化學需氧量的相關系數(shù)均大于0.76，納什效率系數(shù)也在0.58以上，說明在對流擴散模塊中模型參數(shù)設置合理。

表2 污染物變量統(tǒng)計分析

4 情景預測

基于水動力及對流擴散模塊進行驗證后，模型可信度高，該模型能很好模擬出污染物在河道中變化情況。在未來海平面上升或者復雜環(huán)境下，如上游水資源利用量加大或熱帶強風暴潮極端天氣的影響下，對污染物分布進行研究。

通過海南省水功能區(qū)報告，研究區(qū)域設有瓊山工業(yè)、農(nóng)業(yè)用水區(qū)、海口景觀娛樂、漁業(yè)用水區(qū)，水體中總磷含量超出C類景觀娛樂用水水質標準限值(0.05 mg/L)。根據(jù)斷面監(jiān)測顯示，總磷在河口段常年呈持續(xù)超標狀態(tài)，此處選取7月份風暴潮高頻期，上游平均流量207 m3/s，下游平均潮位2.44 m。其他污染分布規(guī)律相似，故以總磷為主要典型污染物進行分布擴散研究。設置以下幾種工況進行對比分析。

a)隨著工農(nóng)業(yè)的發(fā)展，上游需水量不斷加大，考慮上游流量減少20%情況下，量化分析總磷的擴散分布。

b)重現(xiàn)在一次熱帶強風暴潮情況下，即下游潮位上升0.4 m時，量化分析總磷擴散分布。

c)在上游流量減少20%情況下，且考慮下游在強風暴潮極端條件共同作用下，量化分析總磷擴散分布。

5 結果分析與討論

在2016年7月，在3種不同工況下，對主要污染物總磷進行模擬分析，污染物擴散情況見圖8，從圖中可以看出，污染物濃度分布及變化相似，從排污口開始，最大濃度為0.14 mg/L，由于受地形和邊界的影響，污染物主要富集在河道左岸，污染物達到最上游時，濃度約0.02 mg/L。

a)對照組

a)工況一。在上游流量從207 m3/s減少至165.6 m3/s時，下游潮汐作用表現(xiàn)更加明顯，相對對照組來看，橫溝河及海甸溪潮波作用加強，使污染物呈現(xiàn)單條帶狀向上游擴散，由于來流流速變大，污染物濃度迅速降低，但可以看出，污染物擴散距離最遠。

b)工況二?？紤]熱帶風暴潮使下游潮位持續(xù)抬升0.4 m下，上游來流抵達分汊河口處時，強潮流與強來流共同作用下，使得在排污口處，流態(tài)反而變得相對緩慢，流速降低，故污染物主要以擴散機制向上游擴散，擴散距離相對較短，但污染物平均污染帶上濃度較高。

c)工況三。該工況考慮在工況一及工況二2種環(huán)境下污染物擴散問題，模擬結果符合推測，即該結果介于2種工況之間?？梢钥闯?，在工況二的基礎上減少上游來流，潮流作用更加明顯，可以看到污染物以團狀形態(tài)推向上游。

對不同情景下，對污染擴散作定量分析，結果見表3。當上游流量減少20%條件下，污染物將繼續(xù)向上擴散0.6 km，受污染面積減小45%；當下游潮位升高0.4 m時，污染物向上游擴散距離在原始工況下降低22%，污染帶寬度增加60 m，在潮位升高的同時減少上游來流，污染物繼續(xù)向上游擴散。對比工況一和工況二來看，上游流量減少，潮汐作用更加明顯，對于熱帶風暴潮使潮位抬升時，并不能使污染物繼續(xù)向上游擴散，反而會有滯留作用。

表3 不同工況下污染物分布量化統(tǒng)計

6 結論

a)基于非結構網(wǎng)格下連續(xù)性方程和動量方程，建立了南渡江河口水動力及對流擴散模型，通過對流速、流向及潮位進行驗證。通過選取相關系數(shù)及納什效率系數(shù)作為評價指標，結果表明，相關參數(shù)設置可以有效對河道流場進行模擬?；谒畡恿A上，耦合對流擴散模塊，對各污染物降解及擴散系數(shù)進行驗證，通過相關系數(shù)及納什效率系數(shù)分析，表明模型參數(shù)及各污染物系數(shù)選取合理。

b)針對極端天氣條件下最不利熱帶風暴潮的選擇，以模型為基礎，通過升高下游潮位模擬總磷擴散情況，研究表明：升高潮位后，在河口分汊處，上游來流與潮汐作用相互牽制，使得該處流速減小，對流作用較小，污染物主要以擴散形式向上游擴散。

c)在考慮上游用水量加大，模擬總磷擴散情況，結果表明：上游流量減少20%，潮汐作用明顯加強，污染物在潮流推動下繼續(xù)向上游擴散。當同時考慮減少上游流量增加下游潮位時，污染物擴散效果介于兩者單獨作用之間。

d)根據(jù)總磷在不同情景下分布情況，由于在該河口其他污染物屬性與總磷類似，故可參照總磷擴散情況作類推，入?？谥领`山鎮(zhèn)主要是景觀娛樂及漁業(yè)用水區(qū)，水質要求較高，在未來各種環(huán)境下，尤其加強對排污口周圍水環(huán)境監(jiān)測加強排污口污染物排放標準，合理利用上游水資源，優(yōu)化上游用水調度。