潮-徑相互作用下甬江總氮污染的降解、擴散及分布模擬

陳勤思,胡 松,王曉華,劉鵬霞,劉 星,高照詮(.生態(tài)環(huán)境部太湖流域東海海域生態(tài)環(huán)境監(jiān)督管理局生態(tài)環(huán)境監(jiān)測與科學研究中心,上海 005；.上海海洋大學,海洋生態(tài)系統(tǒng)與環(huán)境實驗室,上海 006；.浙江省海洋生態(tài)環(huán)境監(jiān)測中心,浙江 舟山 60)

甬江位于浙江省寧波市,流域面積4518km2,由奉化江和姚江在寧波市區(qū)三江口匯合形成甬江干流,并于鎮(zhèn)海口流入東海.甬江流域是浙江省城市化程度高、經(jīng)濟發(fā)達地區(qū)之一.隨著經(jīng)濟發(fā)展和城市化進程的推進,甬江流域水體污染的問題逐漸顯現(xiàn),臨近海域水質(zhì)污染也十分嚴重.總氮作為甬江入海污染物總量控制的主要指標,其入海行為及影響受到廣泛關(guān)注[1-2].

目前,地表水與海水的監(jiān)測體系相對獨立[3],很難直接依賴監(jiān)測數(shù)據(jù)得到入海總氮污染物的去向和分布情況,數(shù)值模擬成為了解其遷移-轉(zhuǎn)化過程的有效手段.由于不受采樣時間和空間限制,數(shù)值模擬方法可以展示污染物隨水體運動的全過程,常常被用于評估入海污染物對近岸水質(zhì)的影響[4-7].以往研究發(fā)現(xiàn),甬江流域的入海污染物主要輸入到寧波-舟山海域,其擴散和分布情況受徑流[8]、潮汐[9-11]等動力條件影響較大.曹欣中等[12]在對寧波-舟山海域潮流進行模擬的基礎(chǔ)上,分析了污染物隨潮流的輸運特征.張燕等[13]采用河口、陸架和海洋模式(ECOM)對寧波-舟山海域的潮流進行數(shù)值模擬,從而開展了污染物的分布和環(huán)境容量研究.蔡惠文等[14]通過對寧波-舟山海域污染物擴散的數(shù)值模擬,發(fā)現(xiàn)該海域較高的污染物濃度受入海徑流輸入影響較大.這些工作為寧波-舟山海域環(huán)境綜合治理及污染物排海總量控制等提供了支持,具有一定的應(yīng)用價值,但其模擬結(jié)果仍較為理想化,對污染物的分布和影響因素等仍缺少定量的分析.

本文采用數(shù)值模擬方法對甬江入海污染物展開了研究.首先,基于非結(jié)構(gòu)有限體積法海洋模式(FVCOM)建立了寧波-舟山海域的動力場,分析了典型情況下徑流、潮汐對甬江入海污染物分布的影響.進而根據(jù)甬江上游、下游斷面水體總氮濃度的監(jiān)測數(shù)據(jù),確定潮-徑相互作用下甬江總氮污染物的降解系數(shù).最后,在FVCOM 的示蹤物DYE 模塊中增加污染物降解算法,以2021年夏季(7～9月)為例,結(jié)合水體總氮與流量資料,對甬江總氮入海污染物的降解、擴散過程及濃度分布結(jié)果進行了定量分析.本文旨在為該地區(qū)污染物的排放總量控制和陸-海環(huán)境統(tǒng)籌綜合治理服務(wù).

1 數(shù)據(jù)和方法

1.1 動力模型及驗證

采用FVCOM 建立了寧波-舟山海域的水動力模型.FVCOM是基于三維海洋原始方程的非結(jié)構(gòu)網(wǎng)格、有限體積的海洋模式[15],結(jié)合了有限元法擬合邊界、局部加密和有限差分法計算高效等優(yōu)點,能夠并行計算[16],被廣泛應(yīng)用到全球及近海區(qū)域水動力模擬[17].如圖1a所示,本文建立了寧波-舟山海域的數(shù)值模型,模型共生成37959 個三角元與20272 個節(jié)點,其模擬區(qū)域為120.9°～123.3°E、29.25°～31.02°N,包含了杭州灣、寧波-舟山和象山港海域.為保證甬江及寧波-舟山海域網(wǎng)格質(zhì)量,該海域的岸線和地形數(shù)據(jù)采用中華人民共和國海事局出版的海圖予以了訂正,并對該海域進行了網(wǎng)格加密,使海域網(wǎng)格最高水平空間分辨率在30m 以內(nèi).

圖1 寧波-舟山海域FVCOM 模型及驗證Fig.1 FVCOM model grid and tidal calibration in Ningbo-Zhoushan sea area

采用冷啟動方式,模型開邊界由1/30°分辨率的TPXO9 主要分潮的潮汐水位驅(qū)動[18-19],運行時間為2022年3月1日～10月31日,結(jié)果輸出開始時間為2022年3月16日,得到寧波-舟山海域的正壓潮模擬結(jié)果.利用圖1a 中S1 位置2022年4月19日00:00 時至4月21日23:00 時的實測潮位數(shù)據(jù),S2 和S3 位置2022年4月16日07:00 時至4月17日09:00 時的實測潮流數(shù)據(jù)(包含流速、流向)對模型模擬結(jié)果進行了驗證,結(jié)果如圖1b～f所示.將模擬結(jié)果在時間上線性插值到采樣時間后,通過Willmott[20]所提出的方法對模型模擬結(jié)果進行評估,表達式如下:

式中:Pi是模擬值;Oi是實測值;是實測平均值.通過d 值判斷模型效率,d>0.65 時,模型效率極好.經(jīng)驗證,本文模型中潮位、潮流流速及流向的平均d 值分別為0.99,0.91 和0.88,模擬結(jié)果較好.

在部分試驗中,除開邊界的潮位強迫外,模型中還輸入了錢塘江、甬江兩條河流的徑流通量.在動力要素分析部分,徑流流量數(shù)據(jù)參考了《中國河流泥沙公報(2021)》《浙閩臺河流水文資料》[21-22]中發(fā)布錢塘江、甬江的多年平均徑流量.在甬江入海總氮的降解和擴散模擬部分,采用了2021年甬江豐水期(7～9月)平均流量.

1.2 示蹤物模塊及降解算法

非結(jié)構(gòu)有限體積法框架下的示蹤物模塊,起源于全球海洋生態(tài)系統(tǒng)動力學項目(GLOBEC)中對美國喬治淺灘區(qū)域營養(yǎng)鹽輸運的研究,最初由Hu 等[23]建立并用于量化淺灘海域沿鋒面和跨鋒面的深海營養(yǎng)鹽輸運.因設(shè)計初衷是用于驗證比對模型結(jié)果和海上染料投放試驗實測濃度,故取名為DYE 模塊.此后,由于DYE 模塊在物質(zhì)輸運、水體交換等方面具有簡便快速模擬的功能,很快被國內(nèi)外許多研究者應(yīng)用到污染物輸運擴散、水體交換模擬等更多領(lǐng)域的研究中.DYE 模塊利用水動力模塊計算出的三維流速和混合系數(shù)來推算示蹤物濃度的變化.示蹤物被假定為保守量,遵循以下方程:

式中:C 為試驗中保守示蹤物的濃度;D 為總水深;u,v和ω為對應(yīng)x,y 和σ方向流速分量;Kh為垂直擴散系數(shù);FC水平擴散項;CO為保守示蹤物的初始濃度.

在本文的模擬中,DYE 模塊被用于模擬由河流輸入到海水中總氮污染物(營養(yǎng)鹽)的分布.考慮到營養(yǎng)鹽在水體中受物理、化學和生物過程的共同作用下濃度會發(fā)生衰減,并且這類影響對水體中污染物濃度的影響較大[24],因此采用表觀降解系數(shù)方法在DYE 模塊中補充了示蹤物濃度衰減機制,并使其耦合在模型中.表觀降解系數(shù)方法采用降解系數(shù)來表示物質(zhì)衰減的速率,其表達式為:

式中:C0和C 分別表示物質(zhì)降解前后的濃度,mg/L;T為物質(zhì)的降解時間,d;k 即為表觀降解系數(shù),d-1,是在物理、化學和生物等多種因素作用下的綜合降解系數(shù),在本文中表示水體中總氮每日降解的比例.

1.3 總氮數(shù)據(jù)

本文所采用的水體總氮數(shù)據(jù)來自地表水與海水監(jiān)測兩部分.甬江的總氮入海通量根據(jù)歷史流量和總氮數(shù)據(jù)計算得到.在數(shù)值模擬部分采用的地表水監(jiān)測數(shù)據(jù)源于入海河流水質(zhì)監(jiān)測斷面數(shù)據(jù),調(diào)查時間為2021年豐水期(7～9月).寧波-舟山海域海水中總氮數(shù)據(jù)來源于2021年夏季航次近海海水水質(zhì)調(diào)查監(jiān)測項目,調(diào)查時間也為7～9月,以保證數(shù)據(jù)在時間上的匹配性.

2 結(jié)果與討論

2.1 水體總氮的背景調(diào)查

從歷史資料來看(圖2)所示,2020年至2022年間,甬江總氮入海量年際變化非常大.2021年甬江總氮入海量最大,達到2 萬t 以上,并且較前一年(2020年)增加約82%.而在2021年中(圖3),甬江總氮入海量又以第三季度(7～9月)為最高,達到1.21 萬t,約占全年總氮入海量的60%以上.伴隨著這一時期總氮入海量劇增情況的發(fā)生,在同期的海水水質(zhì)監(jiān)測中也觀察到了海水總氮濃度較高的情況.

圖2 甬江總氮入海量年際變化分布Fig.2 Interannual variation of total nitrogen flux in Yongjiang

圖3 2021年甬江總氮入海量季度變化分布Fig.3 Seasonal variation of total nitrogen flux in Yongjiang in 2021

根據(jù)2021年夏季航次近海海水水質(zhì)監(jiān)測數(shù)據(jù),得到了寧波-舟山海域水體中總氮分布結(jié)果.如圖4所示,2021年7～9月寧波-舟山海域水體中總氮濃度整體較高,近岸區(qū)域以0.5mg/L 以上為主.區(qū)域內(nèi)濃度分布整體呈現(xiàn)近岸向近海遞減、北岸向南岸遞減的態(tài)勢.但是,寧波市鎮(zhèn)海區(qū)、北侖區(qū)近岸海域存在一個特別的高值區(qū)域(圖4 方框所示區(qū)域),區(qū)域內(nèi)濃度超過1.0mg/L,海水監(jiān)測點位(P3)最高濃度可達到1.2mg/L 以上,高于其北部海域水體總氮濃度.因此,對圖4 中P1(甬江干流)和P2(甬江入海口)兩處斷面同期的總氮濃度進行分析,發(fā)現(xiàn)同期P1 斷面總氮濃度均值為3.82mg/L,P2 斷面均值為3.01mg/L,甬江輸入到該海域的總氮污染物通量較大,有可能造成區(qū)域性的高值出現(xiàn).

圖4 2021年夏季寧波-舟山海域水體中總氮濃度和監(jiān)測點位分布Fig.4 Total nitrogen concentration and monitoring site distribution in Ningbo-Zhoushan sea waters in summer 2021

2.2 動力要素分析

采用數(shù)值模擬方法進一步評估甬江輸入性污染物對寧波-舟山海域的影響,首先需要分析該海域不同動力因素影響下污染物輸運擴散的特征.由于研究區(qū)域的河流流量在豐水期會有顯著的增加,其增量幅度在每年也會有較大的差異.因此,為了使模擬結(jié)果更具有代表性,可以展現(xiàn)河流在一般通量狀態(tài)下的典型情況,在此部分將輸入河流設(shè)置為甬江、錢塘江的多年平均徑流量.然后,基于FVCOM 模式中經(jīng)典的DYE 模塊,在模型中加入示蹤物,釋放位置對應(yīng)為圖4 中P1 和P2兩處,釋放濃度為1.如表1所示,本文設(shè)計了不同的試驗方案.

表1 數(shù)值模擬試驗方案Table 1 Experimental scenarios of numerical simulation

試驗#1 為僅考慮甬江徑流輸入的情景,用于評估甬江徑流輸入對污染物的擴散和分布的影響.圖5展示了試驗#1條件下釋放示蹤物第60d的分布情況.從模擬結(jié)果中看出,在不考慮潮汐作用的情況下,甬江的輸入使污染物主要向東南方向輸運,最遠可到達大榭島周圍.同時,徑流的沖擊也使污染物存在北向和東向的輸運,但輸運距離較短,且并未到達金塘島沿岸.污染物濃度分布顯示,高值區(qū)域位于甬江河口東南方向并以其為中心向外部遞減.

圖5 試驗#1 情景下示蹤物第60d 分布Fig.5 Distribution of dye in 60d under #1experimental scenarios

僅考慮潮汐作用,圖6b 和c 分別展示了試驗#2 情景下示蹤物第60d 在平潮和停潮時刻污染物的分布情況,圖6a 展示了當日一個潮周期內(nèi)的疊加結(jié)果.與圖5 相比,寧波-舟山海域內(nèi)強潮流的作用使污染物擴散的范圍明顯增大,平潮時污染物北向輸運可達到30.1°N,停潮時污染物東南向輸運可達到122.1°E、29.9°N.疊加結(jié)果顯示,受潮汐作用,污染物影響范圍可到達金塘島南岸.試驗#2 證明了在寧波-舟山海域,潮流是污染物輸運擴散的重要動力因素.

圖6 試驗#2、試驗#3 和試驗#4 情景下示蹤物第60d 分布Fig.6 Distribution of dye in 60d under #2,#3and #4experimental scenarios

在試驗#3 中增加甬江徑流強迫.圖6e 和f 分別展示了試驗#3情景下示蹤物第60d在平潮和停潮時刻污染物的分布情況,圖6d 展示了當日一個潮周期內(nèi)的疊加結(jié)果.與試驗#1 和試驗#2 相比,在潮-徑相互作用下,污染物分布在平潮、停潮和疊加情景下擴散的范圍整體增大.比較圖6a 與圖6d,甬江徑流使污染物影響范圍可到達金塘島西岸,但高值區(qū)域分布變化較小.可見,在該區(qū)域污染物的擴散分布也受到徑流一定程度的影響.

試驗#4 中增加了錢塘江徑流,但其結(jié)果與試驗#3 相比變化不大.試驗#4 證明錢塘江徑流對該區(qū)域污染物輸運和分布影響較小,因此寧波近岸海域總氮高值區(qū)域來源于甬江輸入的可能性較大.

2.3 甬江入海總氮的降解和擴散模擬

在FVCOM 模型中,動力部分采用甬江2021年豐水期(7～9月)平均徑流量和同期開邊界潮位作為強迫場,采用表觀降解系數(shù)方法在FVCOM 模式DYE 模塊中增加衰減機制,以模擬甬江入海總氮污染物的擴散分布情況.根據(jù)圖4 中展示的甬江P1 與P2 斷面位置,在模型中P1 對應(yīng)位置設(shè)置污染物輸入初始濃度為3.82mg/L,根據(jù)P2 位置的污染物濃度3.01mg/L 計算得到式(3)中表觀降解系數(shù)k 約為0.004.

從短期來看,如圖7所示,展示了該情景下第7d總氮污染物的分布情況.圖7b 為第7d 平潮時刻總氮污染物的位置,其顯示漲潮過程中總氮污染物主要影響了鎮(zhèn)海區(qū)東部和北部海水水質(zhì).而在停潮情況下,如圖7c所示,總氮污染物主要分布北侖區(qū)北部海域和甬江口北部一線.圖7a 展示了第7d 一個潮周期內(nèi)總氮污染物分布的疊加結(jié)果,其顯示,在2021年豐水期(7～9月)情景下,甬江總氮污染物輸入到寧波-舟山海域后,短期內(nèi)顯著影響了寧波市鎮(zhèn)海區(qū)和北侖區(qū)近岸海域水質(zhì).從模擬結(jié)果來看,0.05mg/L 包絡(luò)線向北抵達金塘島,向東最大影響到北侖區(qū)穿山半島以東海域,而高值區(qū)主體仍位于寧波市鎮(zhèn)海區(qū)近岸海域.

圖7 2021年7～9月甬江排海總氮污染物第7d 分布Fig.7 Distribution of total nitrogen pollutants discharged from Yongjiang on the 7d from July to September 2021

從長期變化來看,圖8 為2021年夏季(7～9月)甬江排海總氮污染物第15,30,60和120d一個潮周期內(nèi)總氮污染物分布的疊加結(jié)果,其展示了15～90d總氮污染物分布的變化情況.結(jié)果顯示,在2021年豐水期甬江徑流和潮汐的共同作用下,甬江輸入的總氮污染物在15d 內(nèi)影響了寧波-舟山海域中大部分區(qū)域,0.05mg/L 包絡(luò)線向東北方向可到達舟山市本島,向東南方向則擴散到達穿山半島以南海域.30～60d 結(jié)果表明,在近2 個月的時間內(nèi),甬江輸入的總氮污染物可以逐漸影響寧波市和舟山市之間大部分海域,并在60d 使穿山半島周圍海水總氮濃度整體增大約0.1mg/L,使寧波市近岸-金塘島一線和大榭島周圍海域海水總氮濃度整體增大0.2mg/L 以上.60d 與90d 結(jié)果對比表明,考慮總氮的降解過程后,總氮濃度大于0.5mg/L 的高值區(qū)域分布隨時間的變化已較小,而低濃度區(qū)域分布仍有向外延伸的趨勢.據(jù)該情景下第90d(約3 個月)結(jié)果估算,甬江輸入至寧波市鎮(zhèn)海區(qū)東部、北部海域的總氮污染物會使該區(qū)域監(jiān)測站位(圖4 中P3 位置)總氮濃度最高增加0.7mg/L,占該點位監(jiān)測得到總氮濃度值的60%以上.因此,本文認為甬江輸入性的總氮污染來源是2021年夏季該區(qū)域水體總氮濃度出現(xiàn)高值的重要因素.

圖8 2021年7～9月甬江排海總氮污染物15 至90d 分布變化情況Fig.8 Distribution of total nitrogen pollutants discharged from Yongjiang during the 15-90d from July to September 2021

3 結(jié)論

3.1 2021年7～9月,甬江輸入到寧波-舟山海域的總氮污染物可能是造成寧波市鎮(zhèn)海區(qū)、北侖區(qū)近岸海域海水總氮濃度出現(xiàn)高值的重要因素.

3.2 潮汐和甬江徑流都會影響寧波-舟山海域污染物的擴散輸運過程.在同時考慮潮-徑相互作用的情景下,甬江輸入到寧波-舟山海域的總氮污染物將在短期內(nèi)顯著影響到寧波市鎮(zhèn)海區(qū)東部、北部海域海水質(zhì)量.

3.3 根據(jù)2021年7～9月的模擬結(jié)果,考慮降解過程后,總氮污染物在15d 內(nèi)影響了寧波-舟山海域中大部分區(qū)域,總氮濃度大于0.5mg/L的高值區(qū)域分布在60～90d 之間變化已較小.根據(jù)第90d 結(jié)果估算,甬江輸入至寧波市鎮(zhèn)海區(qū)東部、北部海域的總氮污染物會使該區(qū)域海水監(jiān)測站位總氮濃度最高增加0.7mg/L.