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

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

(1.西華大學(xué)能源與動(dòng)力工程學(xué)院，四川 成都 610039；2.西華大學(xué)流體及動(dòng)力機(jī)械教育部重點(diǎn)實(shí)驗(yàn)室，四川 成都 610039)

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

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

根據(jù)2018年海南省海洋環(huán)境狀況公報(bào)可知，在2018年有9個(gè)臺(tái)風(fēng)進(jìn)入南海或在南海生成，其中，1804號(hào)臺(tái)風(fēng)(熱帶風(fēng)暴級(jí))“艾云尼”和1809號(hào)臺(tái)風(fēng)(熱帶風(fēng)暴級(jí))“山神”登陸海南島，在海南島沿岸引發(fā)不同程度的風(fēng)暴增水，根據(jù)潮位站資料顯示，海口秀英驗(yàn)潮站最大增水0.38 m，瓊海博鰲驗(yàn)潮站最大增水0.64 m。根據(jù)《海南省水資源公報(bào)》可知，海口市用水量加大，主要是海口常駐人口逐年遞增，工農(nóng)業(yè)用水量加大，上游需水量逐年遞增。隨著海口市城市化進(jìn)程的不斷加快，河道面積縮減、城市排水不暢、水質(zhì)惡化等水環(huán)境問題日趨嚴(yán)重[9]。

為此，根據(jù)南渡江河口水環(huán)境情況，對(duì)南渡江河口水質(zhì)變化規(guī)律進(jìn)行研究。基于潮汐和上游來流相互作用下，開展在不同環(huán)境下污染物在河道中變化規(guī)律。實(shí)現(xiàn)在熱帶復(fù)雜河口基礎(chǔ)上，探究極端天氣和惡劣條件下污染物變化情況，提出水生態(tài)保護(hù)與綜合治理建議和方案，為海南省“十四五”規(guī)劃水環(huán)境治理提供依據(jù)。

1 研究區(qū)域概況

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

南渡江入海口常年平均水體溫度在25℃，有著獨(dú)特的臺(tái)風(fēng)季節(jié)特征，南渡江流域降雨量分配非常不均，每年5—10月降雨量就達(dá)到全年80%，河道水量主要靠降雨補(bǔ)給。

2 研究方法

2.1 MIEK21模型基本原理

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

圖1 南渡江河口研究范圍

2.1.1控制方程

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

(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——風(fēng)速及在x、y方向上的分量，m/s；f——風(fēng)阻力系數(shù)；τxx、τxy、τyy——剪切應(yīng)力分量，kg/m2；pa——當(dāng)?shù)卮髿鈮簭?qiáng)，Pa。

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

(4)

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

3 模型構(gòu)建

3.1 網(wǎng)格剖分

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

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

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

3.2 地形插值

利用項(xiàng)目區(qū)河口段1∶100 00實(shí)測(cè)水下地形圖，提取相應(yīng)河道高程點(diǎn)，基于非結(jié)構(gòu)性網(wǎng)格基礎(chǔ)上，導(dǎo)入河道高程文件，采用自然鄰近插值法進(jìn)行，得到河底地形(圖3)。

圖3 河床高程分布(m)

3.3 邊界條件及源匯設(shè)置

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

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

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

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

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

3.4 重要參數(shù)設(shè)置

3.4.1糙率

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

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

圖5 糙率場(chǎng)分布(m)

3.4.2擴(kuò)散系數(shù)

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

3.4.3降解系數(shù)

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

Ct=Coe-kt

(5)

(6)

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

根據(jù)該公式初估各種污染物的降解系數(shù)，再經(jīng)過反復(fù)驗(yàn)證后，最終確定南渡江河口區(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 模擬驗(yàn)證

模型水動(dòng)力計(jì)算時(shí)間選取在夏季豐水時(shí)期，時(shí)間從2016年6月28日19:00至2016年6月29日23:00。通過入海口分汊點(diǎn)L1實(shí)測(cè)值，對(duì)模型的流速、流向和水位進(jìn)行驗(yàn)證。流速、水位及流向(驗(yàn)證)見圖6。

a)流速

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

(7)

(8)

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

3.6 MIKE 21 AD模擬及驗(yàn)證

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

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

4 情景預(yù)測(cè)

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

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

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

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

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

5 結(jié)果分析與討論

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

a)對(duì)照組

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

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

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

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

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

6 結(jié)論

a)基于非結(jié)構(gòu)網(wǎng)格下連續(xù)性方程和動(dòng)量方程，建立了南渡江河口水動(dòng)力及對(duì)流擴(kuò)散模型，通過對(duì)流速、流向及潮位進(jìn)行驗(yàn)證。通過選取相關(guān)系數(shù)及納什效率系數(shù)作為評(píng)價(jià)指標(biāo)，結(jié)果表明，相關(guān)參數(shù)設(shè)置可以有效對(duì)河道流場(chǎng)進(jìn)行模擬。基于水動(dòng)力基礎(chǔ)上，耦合對(duì)流擴(kuò)散模塊，對(duì)各污染物降解及擴(kuò)散系數(shù)進(jìn)行驗(yàn)證，通過相關(guān)系數(shù)及納什效率系數(shù)分析，表明模型參數(shù)及各污染物系數(shù)選取合理。

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

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

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