基于水動(dòng)力數(shù)值模擬的羅源灣水體交換周期及納潮量計(jì)算

陳志琦，逄 勇，張 倩，孫 凡，吳昌淦

(1.河海大學(xué)淺水湖泊綜合治理與資源開發(fā)教育部重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210098；2.河海大學(xué)環(huán)境學(xué)院，江蘇 南京 210098)

隨著經(jīng)濟(jì)的飛速發(fā)展，環(huán)境污染問(wèn)題日益突出[1-2]。沿海城市在依靠?jī)?yōu)越的地理位置發(fā)展沿海經(jīng)濟(jì)的同時(shí)，廢水和廢熱不斷排入海灣，影響了海灣的生態(tài)環(huán)境。羅源灣位于福州市東北部濱岸處，隨著羅源灣灘涂圍墾面積的增加和臨港工業(yè)及養(yǎng)殖業(yè)發(fā)展，羅源灣水質(zhì)逐步惡化，生態(tài)系統(tǒng)的敏感性越來(lái)越強(qiáng)。近年來(lái)，福建省政府對(duì)羅源灣開展了退養(yǎng)工作，羅源灣水質(zhì)雖有所改善[3-4]，但仍存在超標(biāo)現(xiàn)象。水體交換周期及納潮量是衡量海灣自凈能力及海灣生存能力的重要指標(biāo)[5-6]，研究海灣的納潮量及水體交換周期對(duì)海域的合理開發(fā)利用具有重要意義。

目前，對(duì)于水體交換周期的研究主要采用箱式模型、粒子追蹤模型和對(duì)流-擴(kuò)散水質(zhì)模型。蔣磊明等[7]根據(jù)FVCOM(finite-volume coastal ocean model)模擬結(jié)果計(jì)算了欽州灣的納潮量及水體交換周期。秦曉等[8]基于Mike21軟件對(duì)東山灣進(jìn)行了水動(dòng)力數(shù)值模擬，并計(jì)算得到東山灣納潮量及水體交換周期。周爭(zhēng)橋等[9]采用Lagrange質(zhì)點(diǎn)追蹤方法對(duì)防城灣豐、枯水期水體交換能力進(jìn)行了評(píng)價(jià)。1988年開始，國(guó)內(nèi)學(xué)者對(duì)羅源灣開展了研究，但多集中在灣內(nèi)魚貝類等水生生物的研究，對(duì)羅源灣水動(dòng)力模擬較少。胡建宇[10]根據(jù)實(shí)測(cè)資料及單箱模型計(jì)算公式計(jì)算得出羅源灣納潮量較大，大潮期間納潮量可達(dá)9.6億m3，并假設(shè)羅源灣灣外海水與灣內(nèi)海水均勻混合，得出羅源灣水體半交換周期約為17個(gè)潮周期。杜伊等[11-12]應(yīng)用ECOMSED(estuarine, coastal and ocean modeling system with sediments)模擬了羅源灣水動(dòng)力情況，并用示蹤粒子三維追蹤模塊分析羅源灣水體交換情況，結(jié)果表明，羅源灣內(nèi)不同海區(qū)海水存留時(shí)間存在差異，主要由歐亞余流場(chǎng)結(jié)構(gòu)差異所致；西南風(fēng)有利于水體交換，東北風(fēng)阻礙水體交換，高潮及落潮期間水體交換能力強(qiáng)，表層水體交換快，底層較慢。王勇智等[13]采用Mike3軟件，分析了羅源灣3個(gè)典型圍填海時(shí)期岸線與地形條件下的納潮量和水體交換能力變化，結(jié)果表明，1996年和2012年的海灣全潮平均納潮量與20世紀(jì)60年代相比分別減少了約20.59%和28.38%，灣內(nèi)30 d的平均水體交換率則分別減小了19.17%和21.42%，水體半交換時(shí)間延長(zhǎng)了約1.74 d和2.42 d。

納潮量是海灣環(huán)境評(píng)價(jià)的重要指標(biāo)，近幾年，國(guó)內(nèi)學(xué)者應(yīng)用Mike等數(shù)值模擬軟件對(duì)國(guó)內(nèi)外各個(gè)海域的納潮量進(jìn)行了模擬分析。婁安剛等[14]運(yùn)用二維潮波運(yùn)動(dòng)方程建立了丁字灣附近海域的水動(dòng)力模型，模擬規(guī)劃實(shí)施前后水動(dòng)力特性及納潮量的變化；陳靜等[6]應(yīng)用Mike21軟件計(jì)算了大連灣及3個(gè)內(nèi)灣的動(dòng)態(tài)納潮量；王誠(chéng)超等[5]基于樂(lè)清灣近50年水域面積的演變特征，對(duì)樂(lè)清灣納潮量的演變規(guī)律進(jìn)行了研究，得出樂(lè)清灣納潮量存在遞減趨勢(shì)，納潮量對(duì)水域面積改變的響應(yīng)關(guān)系有階段性特征。

目前，應(yīng)用數(shù)值模擬軟件模擬羅源灣水動(dòng)力情況，并計(jì)算其水體交換周期以及納潮量的研究較少，本文應(yīng)用Mike21軟件模擬羅源灣水動(dòng)力特性，計(jì)算和分析羅源灣納潮量及水體交換周期，可為羅源灣污染物總量控制提供基礎(chǔ)性資料。

1 研究區(qū)概況

羅源灣是一個(gè)半封閉港灣，口窄腹大。海灣面積約154 km2，位于東經(jīng)119.60°～119.86°、北緯26.33°～26.49°之間，北部與羅源縣相連，南部與連江縣毗鄰，羅源縣縣城與羅源灣西北角相接。匯入羅源灣的地表徑流主要有羅源縣的起步溪、護(hù)國(guó)溪、南門溪，連江縣的鯉溪、牛柘溪等，同時(shí)，羅源灣沿海陸域因受地形自然切割的影響，形成眾多短促的獨(dú)自入海溪流。灣內(nèi)有1個(gè)潮位站——跡頭站，羅源灣地理位置及地形如圖1所示。

圖1 羅源灣地理位置及地形

2 研究方法

主要采用Mike21軟件的水動(dòng)力模塊和對(duì)流擴(kuò)散模塊進(jìn)行數(shù)值模擬。Mike21軟件的計(jì)算基于三向不可壓縮和Reynolds值均勻分布的Navier-Stokes方程，并采用靜水壓力的假定和Boussinesq假定[15]。其水動(dòng)力模塊基本計(jì)算方程包括連續(xù)性方程和動(dòng)量方程，對(duì)流擴(kuò)散模塊采用二維對(duì)流擴(kuò)散方程進(jìn)行計(jì)算。本文根據(jù)模型計(jì)算的潮位值確定羅源灣納潮量，依據(jù)不同時(shí)刻保守物質(zhì)質(zhì)量濃度確定羅源灣水體交換周期。

2.1 納潮量計(jì)算方法

在一個(gè)潮周期內(nèi)海灣所能接納的海水體積即為海灣的納潮量。為獲得精確的計(jì)算結(jié)果，在有限體積方法中，將納潮量概念應(yīng)用到每個(gè)網(wǎng)格范圍內(nèi)：

(1)

式中:P為納潮量；Si為第i個(gè)網(wǎng)格的面積；Hhi和Hli分別為第i個(gè)網(wǎng)格的高、低潮水位；n為所選取海域的網(wǎng)格總數(shù)。

2.2 水體交換周期計(jì)算方法

設(shè)羅源灣保守物質(zhì)的初始質(zhì)量濃度為1 mg/L，外灣質(zhì)量濃度為0 mg/L，無(wú)其他源項(xiàng)匯入，保守物質(zhì)衰減系數(shù)取0。灣內(nèi)污染物通過(guò)水體交換向?yàn)惩廨斠疲?dāng)灣內(nèi)保守物質(zhì)質(zhì)量濃度降至0.5 mg/L時(shí)，認(rèn)為水體完成了半交換，此時(shí)所用的時(shí)間即為該海灣水體的半交換周期。水體交換率指某點(diǎn)保守物質(zhì)質(zhì)量濃度的減小值與初始質(zhì)量濃度的比值：

(2)

式中:γ為交換率；Δρ為保守物質(zhì)質(zhì)量濃度的減小值，mg/L；ρ1為保守物質(zhì)初始質(zhì)量濃度，mg/L。

3 模型構(gòu)建與參數(shù)率定

模型計(jì)算范圍是羅源灣灣口至羅源灣整個(gè)海域，采用三角形網(wǎng)格對(duì)羅源灣海域進(jìn)行剖分，網(wǎng)格總數(shù)6 821，最小的三角形單元面積為2 321 m2。模型計(jì)算的時(shí)間步長(zhǎng)為Δt= 300 s，計(jì)算時(shí)間為2019年9月1日至2019年12月31日，總時(shí)長(zhǎng)為122 d。

初始水位取模型計(jì)算起始時(shí)刻邊界潮位均值，設(shè)為-1.7 m，起始時(shí)刻流速設(shè)為0 m/s，降雨、蒸發(fā)、風(fēng)速和風(fēng)向資料來(lái)源于國(guó)家氣象信息中心氣象監(jiān)測(cè)站——寧德站(區(qū)站號(hào)58846)2019年逐日實(shí)測(cè)數(shù)據(jù)。模型開邊界采用潮位驅(qū)動(dòng)，灣口處潮位由Mike21軟件Toolbox中的全球潮汐預(yù)測(cè)程序給出，該潮汐預(yù)測(cè)程序采用預(yù)測(cè)區(qū)域8個(gè)分潮的調(diào)和常數(shù)進(jìn)行疊加獲得潮位，再根據(jù)計(jì)算域內(nèi)潮位驗(yàn)證情況進(jìn)行調(diào)整。岸灘邊界的變化由計(jì)算格點(diǎn)水深來(lái)判別，當(dāng)某個(gè)網(wǎng)格點(diǎn)的水深小于某一極限值(取0.001 m)，且其周圍網(wǎng)格點(diǎn)的水位均低于該網(wǎng)格點(diǎn)的水位時(shí)，該網(wǎng)格點(diǎn)被認(rèn)為是“干單元”，不參與計(jì)算；當(dāng)“干單元”周圍任意一點(diǎn)的水位高于該點(diǎn)水位時(shí)，“干單元”即轉(zhuǎn)化為“濕單元”，重新參與計(jì)算，模型計(jì)算區(qū)域網(wǎng)格劃分如圖2所示。

選取羅源灣跡頭站大小潮潮位對(duì)模型進(jìn)行率定。小潮潮位率定采用2019年9月7—9日潮位資料，大潮潮位率定采用2019年10月27—29日潮位資料。模型對(duì)大潮潮位模擬的平均絕對(duì)誤差為26 cm，平均相對(duì)誤差為14%；模型對(duì)小潮潮位模擬的平均絕對(duì)誤差為14 cm,平均相對(duì)誤差為19%。模型率定結(jié)果如圖3所示。

圖2 羅源灣模型計(jì)算區(qū)域與網(wǎng)格劃分

(a)大潮

(b)小潮

余東等[16]研究表明，羅源灣可門水道處流速約為1.0 m/s，灣內(nèi)流速約為0.5 m/s；漲潮時(shí)，灣外水流經(jīng)灣口可門水道流向?yàn)硟?nèi)后分為兩支，一支流向西北灣跡頭，另一支流向西南；落潮時(shí)，羅源灣海域流場(chǎng)是典型的往復(fù)流，與本文流場(chǎng)模擬結(jié)果基本相符。

根據(jù)率定結(jié)果得到模型灣內(nèi)海域Smagorinsky系數(shù)為0.28，糙率為0.024～0.035，風(fēng)拖曳系數(shù)為0.001 0～0.001 5。

4 模擬結(jié)果與分析

4.1 流場(chǎng)模擬

羅源灣海域流場(chǎng)變化模擬結(jié)果如圖4所示。

a.漲潮時(shí)，在灣外海域潮波作用下，大量海水涌入灣內(nèi)，潮流在可門水道轉(zhuǎn)向西南進(jìn)入羅源灣，灣口處流速較大，最大處可達(dá)1.24 m/s。水流進(jìn)入灣口后，水域面積增大，流速降低，為0.60～1.24 m/s，羅源灣內(nèi)部流速為0.20～0.60 m/s。漲潮時(shí)羅源縣城區(qū)附近海域以及將軍帽工業(yè)園區(qū)附近海域流速受潮波影響較小，低于0.20 m/s。總體來(lái)說(shuō)，羅源灣灣口處水動(dòng)力條件較好，北岸水動(dòng)力條件較差，北岸流場(chǎng)受潮波影響程度小于南岸流場(chǎng)。

(a)大潮(落潮) (b)大潮(漲潮)

(c)小潮(落潮) (d)小潮(漲潮)

b.落潮時(shí)，潮流從羅源灣退至外海，流向與漲潮時(shí)基本相反，海水由北岸流經(jīng)南岸，再?gòu)哪习读飨驗(yàn)晨?。灣?nèi)流速為0.26～0.90 m/s，可門水道內(nèi)流速較大，最大可達(dá)1.56 m/s，水動(dòng)力條件較好，北岸水動(dòng)力條件劣于南岸。落潮時(shí)南北岸流場(chǎng)密集，大量海水從可門水道流出羅源灣。

4.2 水體交換周期

按照水體交換周期研究方法進(jìn)行參數(shù)設(shè)置，計(jì)算可得15 d、30 d、45 d、60 d后羅源灣水體平均交換率分別為48.5%、69.3%、78.8%和85.5%，半交換周期約為15～16 d (表1)。由圖5可知，在相同時(shí)間內(nèi)灣口及可門水道內(nèi)水體交換率最大；灣中部水體交換率中等；灣頂處(羅源灣西北角)水體交換率最小。15 d后灣中部水體交換率在30%～80%之間,30 d后可達(dá)50%～90%;15 d后灣頂處水體交換率低于20%,30 d后交換率可達(dá)44%,60 d 后可達(dá)72%。這表明隨著時(shí)間的推移，受灣外海域潮波影響，灣內(nèi)的保守物質(zhì)濃度不斷降低，其中距離外海最近的灣口及可門水道水體交換能力最強(qiáng)，污染物濃度迅速降低。灣中部水體交換能力良好，而西北角灣頂處水體交換能力最弱，這主要是由于羅源灣西北角三面環(huán)陸,而整個(gè)海灣僅通過(guò)東南方可門水道與外界進(jìn)行水體交換,導(dǎo)致此處受外海潮波影響較小，且流速較低，不利于污染物擴(kuò)散。根據(jù)1998年胡建宇[10]對(duì)羅源灣水體半交換周期的研究，羅源灣水體半交換周期為17個(gè)潮周期，約為9 d；本文得到羅源灣水體半交換周期約為15 d，說(shuō)明羅源灣海域更新能力減弱，半交換周期時(shí)間延長(zhǎng)。

表1 羅源灣水體交換率

(a)15 d(落潮) (b)15 d(漲潮)

(c)30 d(落潮) (d)30 d(漲潮)

(e)45 d(落潮) (f)45 d(漲潮)

(g)60 d(落潮) (h)60 d(漲潮)

為定量分析灣內(nèi)污染物交換情況，在灣口、灣中部、灣西部各取一點(diǎn)位(圖1)，分析其保守物質(zhì)質(zhì)量濃度的變化規(guī)律。由圖6可知，受漲落潮影響，A、B、C點(diǎn)保守物質(zhì)質(zhì)量濃度呈波動(dòng)變化，灣口保守物質(zhì)質(zhì)量濃度受漲落潮影響最為明顯。位于灣口A點(diǎn)的保守物質(zhì)質(zhì)量濃度急劇降低后趨于穩(wěn)定，3 d后A點(diǎn)的保守物質(zhì)質(zhì)量濃度低于0.5 mg/L。位于灣中部B點(diǎn)的保守物質(zhì)質(zhì)量濃度降低速率無(wú)明顯的突變過(guò)程，26 d后B點(diǎn)的保守物質(zhì)質(zhì)量濃度基本低于0.5 mg/L，50 d后B點(diǎn)的保守物質(zhì)質(zhì)量濃度低于0.2 mg/L。位于西北角C點(diǎn)的保守物質(zhì)質(zhì)量濃度的變化受漲落潮影響最小，C點(diǎn)的水體半交換周期約為43 d，60 d后C點(diǎn)的保守物質(zhì)質(zhì)量濃度約為0.3 mg/L，水體交換率僅為70%。

圖6 A、B、C點(diǎn)保守物質(zhì)質(zhì)量濃度隨時(shí)間變化曲線

4.3 納潮量

利用羅源灣模型模擬結(jié)果計(jì)算9—10月(秋季)的納潮量，羅源灣平均納潮量為5.73 億m3，大潮納潮量為7.07億m3，小潮納潮量為4.40億m3；大小潮時(shí)期納潮量差距較大，兩者比值為1.61∶1。根據(jù)胡建宇[10]的研究成果，1998年羅源灣平均納潮量為7.2億m3，大潮納潮量為9.6億m3，小潮納潮量為5.1億m3。與1998年相比，羅源灣納潮量有減少的趨勢(shì)，2019年平均納潮量減少了20.4%，大潮納潮量減少了26.4%，小潮納潮量減少了13.7%。根據(jù)杜伊等[11-12]對(duì)羅源灣納潮量的研究，2007年羅源灣平均納潮量為7.2億m3，與該結(jié)果相比，羅源灣納潮量仍呈減少的趨勢(shì)。

綜合考慮羅源灣海域形狀及水動(dòng)力特點(diǎn)，將羅源灣分為Ⅰ區(qū)、Ⅱ區(qū)和Ⅲ區(qū)，分別計(jì)算其納潮量。Ⅰ區(qū)、Ⅱ區(qū)和Ⅲ區(qū)面積占羅源灣總面積的比例分別為16%、66%和18%。Ⅲ區(qū)位于羅源灣西部，受外海水動(dòng)力條件影響較小，水體交換率較小，15 d后水體交換率低于20%，大潮潮差為5.55 m，小潮潮差為3.42 m，平均納潮量為0.95 億m3；Ⅱ區(qū)范圍較大，水體交換能力中等，大潮潮差為5.15 m，小潮潮差為3.28 m，平均納潮量為3.65 億m3；Ⅰ區(qū)包括灣口及可門水道，受外海區(qū)域潮波影響大，潮差較大，大潮潮差為6.34 m,小潮潮差為3.7 m，平均納潮量為1.13億m3。羅源灣秋季納潮量計(jì)算結(jié)果見表2。

表2 羅源灣秋季納潮量計(jì)算結(jié)果

在模型計(jì)算基礎(chǔ)上，分析秋季(9—10月)、冬季(12月)海灣納潮量變化。結(jié)果表明，冬季Ⅲ區(qū)、Ⅱ區(qū)和Ⅰ區(qū)的平均納潮量分別為秋季的87%、89%和86%，整個(gè)海灣冬季納潮量減少，為秋季的88%。

5 結(jié) 論

a.構(gòu)建了羅源灣水動(dòng)力數(shù)學(xué)模型，率定得到模型灣內(nèi)海域Smagorinsky系數(shù)為0.28，糙率為0.024～0.035，風(fēng)拖曳系數(shù)為0.001 0～0.001 5。

b.漲潮時(shí)，灣口處流速達(dá)到1.24 m/s，海水流過(guò)灣口后，流速降低為0.60～1.24 m/s，羅源灣內(nèi)部流速為0.20～0.60 m/s；落潮時(shí)，海水流向與漲潮時(shí)基本相反，灣內(nèi)流速為0.26～0.90 m/s，灣口處最大可達(dá)1.56 m/s。羅源灣灣口處水動(dòng)力條件較好，北岸水動(dòng)力條件較差，北岸流場(chǎng)受外海潮波影響程度小于南岸流場(chǎng)。

c.羅源灣水體平均半交換周期約15 d，15 d、30 d、45 d和60 d后羅源灣水體平均交換率分別為48.5%、69.3%、78.8%和85.5%，灣口及可門水道處水體交換率最大，交換能力最強(qiáng)，羅源灣西北角和西南處的灣頂處水體交換率最小。

d.秋季時(shí)，羅源灣平均納潮量為5.73億m3，大潮納潮量為7.07億m3，小潮時(shí)納潮量為4.40億m3。Ⅰ區(qū)平均納潮量為1.13億m3，占總納潮量的19.7%；Ⅱ區(qū)平均納潮量為3.65億m3，占總納潮量的63.7%；Ⅲ區(qū)平均納潮量為0.95億m3，占總納潮量的16.6%。

e.羅源灣西北角流速小，水體交換率低，納潮量有限，水動(dòng)力條件較差。建議在羅源灣入海污染物通量管控過(guò)程中，綜合考慮海域自身水動(dòng)力條件，合理控制羅源灣西北角入海污染物排放量。