林偉波，魏愛泓，冒士鳳

(江蘇省海涂研究中心，南京 210036)

灌河是蘇北唯一一條在干流上沒有建閘控制的天然大型河道，支流眾多, 其下游在堆溝至燕尾港處有新沂河匯入;干流西起灌南縣境內(nèi)的東三岔, 東至燕尾港的灌河口, 全長(zhǎng)74.5 km，流域面積6 400 km2[1,2]。灌河潮位、水流主要受黃海潮波控制, 受徑流影響較小。灌河口海區(qū)潮汐屬于非正規(guī)半日潮，潮汐特性表現(xiàn)為潮差大，流急，灌河口內(nèi)段， 因受邊界約束潮波變形， 形成前進(jìn)駐波混合型潮波， 表現(xiàn)為前坡陡、后坡緩，通常漲潮歷時(shí)小于落潮歷時(shí)[3]。研究該區(qū)域的水動(dòng)力特性，對(duì)于控制灌河口海域的污染具有極其重要的意義。Chen等人[4]成功地建立了三維非結(jié)構(gòu)、原始方程、有限體積的海洋模型(FVCOM)，并對(duì)我國(guó)東部海域[5]和美國(guó)部分海灣[6-9]進(jìn)行了大量的數(shù)值計(jì)算。這個(gè)模型在理論和數(shù)值計(jì)算方面基本解決了淺海陸架、河口物理海洋和生態(tài)動(dòng)力學(xué)模型中復(fù)雜岸界擬合和計(jì)算速度的難題，該模型集水動(dòng)力模塊、泥沙輸運(yùn)模塊、生態(tài)模塊和水質(zhì)預(yù)測(cè)模塊一體，可用于模擬水系統(tǒng)二維和三維流場(chǎng)、物質(zhì)輸運(yùn)(包括溫、鹽、非黏性和黏性泥沙的輸運(yùn))、生態(tài)過(guò)程及淡水入流。其模擬范圍包括河口、河流、湖泊、水庫(kù)、濕地以及自近岸到陸架的海域。Lucy 等人[10]基于FVCOM模型研究外海潮流入侵感潮河道受地形、河道寬度、流量和潮汐動(dòng)力相互影響。Mochael Togneri 等人[11]基于ROMS模型研究潮汐紊動(dòng)能量，通過(guò)與ADCP實(shí)測(cè)數(shù)據(jù)比對(duì)，表明ROMS模型在預(yù)測(cè)高能量區(qū)的潮汐動(dòng)力紊動(dòng)能量有很好的效果。Mohammad NabiAllahdad等人[12]基于三維斜壓模型FVCOM研究分層流對(duì)路易安娜陸架潮流動(dòng)力精確模擬的影響及其解決方法。 近幾年FVCOM模式在國(guó)內(nèi)也開始得到應(yīng)用。。宋德海[13]等基于采用不規(guī)則三角網(wǎng)格和有限體積方法的FVCOM 模式，建立欽州灣三維潮流數(shù)值模型來(lái)重現(xiàn)欽州灣的潮位和潮流變化狀況。李希彬[14]等基于FVCOM模型在湛江附近海域建立了三維動(dòng)邊界水動(dòng)力模型，模擬計(jì)算了湛江東海島填海大堤現(xiàn)狀以及1958年大堤修建之前湛江海域的水動(dòng)力場(chǎng)。Xuan Jiliang等人[15]基于FVCOM模型研究長(zhǎng)江口余流組成機(jī)制及其在平均流中的重要作用。本文基于該模型建立了灌河口三維潮流數(shù)學(xué)模型，對(duì)灌河口海域水動(dòng)力過(guò)程進(jìn)行了數(shù)值模擬，模擬結(jié)果和實(shí)測(cè)資料比較吻合，并且對(duì)模擬的流場(chǎng)分布進(jìn)行了比較詳細(xì)的分析，表明該模式可以用于模擬和分析河口以及海洋動(dòng)力場(chǎng)的分布以及變化特征，也為研究灌河口的污染物擴(kuò)散模型提供正確的水動(dòng)力條件，從而保證了預(yù)測(cè)結(jié)果的可靠性。

1 三維有限體積模型及計(jì)算方法

1.1 控制方程

在流體不可壓縮、Boussinesq和靜力近似下,給出雷諾平均的三維紊流河口海岸海洋控制方程組。海洋控制方程組是由動(dòng)量方程、連續(xù)方程、溫度方程、鹽度方程、密度方程和紊流方程組成。

動(dòng)量方程：

(1)

(2)

(3)

連續(xù)方程：

(4)

溫度和鹽度方程：

(5)

(6)

密度方程：

ρ=ρ(T,S)

(7)

式中：x、y、z分別為Cartesian直角右手坐標(biāo)系的東、北和垂直方向坐標(biāo)；u、v、w分別為x、y、z方向上的速度分量；T為海水位溫；S為海水鹽度；P為壓強(qiáng)；f為科式力參數(shù)；g為重力加速度；ρ為海水密度；Km為垂直渦黏性系數(shù)；Kh為垂直熱力擴(kuò)散系數(shù)；Fu、Fv、FT和Fs分別為水平動(dòng)量、熱力和鹽度擴(kuò)散項(xiàng)。

Km和Kh由修正的Mellor和Yamada的2.5階湍流閉合子模型計(jì)算。

(8)

(9)

1.2 邊界條件

u，v，w的一般性表面和底部邊界條件如下。

在表面z=ζ(x,y,t)處：

(10)

而在底部z=-H(x,y)處：

(11)

阻力系數(shù)Cd由下面對(duì)數(shù)底邊界層計(jì)算值和常數(shù)值中的最大值確定，即：

Cd=max[k2/ln(zab/z0)2, 0.002 5]

(12)

式中：Z0為底部粗糙度參數(shù)。

溫度的表面和底部邊界條件如下。

在表面z=ζ(x,y,t)處：

(13)

在z=-H(x,y)處：

(14)

式中：Qn(x,y,t)為表面凈熱通量，它等于短波輻射、長(zhǎng)波輻射、感熱通量和潛熱通量之和；SW(x,y,0,t)為在海表面處短波輻射通量；cp為海水比熱系數(shù)。

鹽度在表面和底部的邊界條件如下。

在表面z=ζ(x,y,t)處：

(15)

在底部z=-H(x,y)處：

(16)

通過(guò)引入底部邊界條件可在模型中考慮地下水通量的作用。湍流動(dòng)能和混合長(zhǎng)度方程組的表面和底部邊界如下。

在表面z=ζ(x,y,t)處：

(17)

在底部z=-H(x,y)處：

(18)

式中：uτs和uτb分別為表面和底部對(duì)數(shù)邊界層的摩擦速度。

側(cè)邊界處運(yùn)動(dòng)學(xué)、熱量和鹽度的邊界條件分別為：

(19)

式中：n為邊界的法向坐標(biāo)軸；vn則為邊界法向速度分量。

1.3 數(shù)值方法

模型采用內(nèi)外模分離法求解。二維外模數(shù)值格式為基于三角形網(wǎng)格的有限體積法，將連續(xù)方程、動(dòng)量方程在三角形區(qū)域積分后，通過(guò)修正過(guò)的四階龍格庫(kù)塔法求解。三維內(nèi)模的動(dòng)量方程的求解采用簡(jiǎn)單的顯式和隱式相結(jié)合的差分格式求解，其中流速的局部變化項(xiàng)采用二階精度的龍格庫(kù)塔時(shí)間積分格式，對(duì)流項(xiàng)采用二階精度的迎風(fēng)格式,垂直擴(kuò)散項(xiàng)則用顯示格式求解。

2 計(jì)算結(jié)果與分析

2.1 驗(yàn)證資料說(shuō)明

驗(yàn)證所用資料為2012年由長(zhǎng)江下游水文水資源勘測(cè)局開展了研究海域的水文調(diào)查，布置3個(gè)潮位觀測(cè)站，在全潮水文測(cè)驗(yàn)前一天開始潮位觀測(cè)，連續(xù)觀測(cè)10 d；布置6條垂線，進(jìn)行大、小潮連續(xù)28 h以上的全潮水文測(cè)驗(yàn)，測(cè)站布設(shè)如圖1所示。測(cè)驗(yàn)內(nèi)容包括10月8日-10月19日的連續(xù)10 d潮位以及一個(gè)連續(xù)潮汛大、小潮潮流，時(shí)間間隔為1 h；小潮(10月8日-10月9日)和大潮(10月18日-10月19日)，28 h船測(cè)垂線流速，每條測(cè)驗(yàn)垂線自海底至水面等間距布設(shè)6點(diǎn)。

圖1 研究海域水文測(cè)站布置Fig.1 The location of hydrometric stations

2.2 模型范圍和參數(shù)設(shè)置

模型的計(jì)算域范圍為119°12′E～121°12′E，33°45′N～35°36′N，近岸海域網(wǎng)格步長(zhǎng)500 m，外海開邊界網(wǎng)格步長(zhǎng)3 km。計(jì)算區(qū)域水深特征如圖2所示,計(jì)算網(wǎng)格如圖3所示。計(jì)算區(qū)域總共由12 931個(gè)網(wǎng)格單元，6 682個(gè)網(wǎng)格節(jié)點(diǎn)組成。水動(dòng)力模型的參數(shù)校核包括底部摩擦系數(shù)、垂直渦黏系數(shù)背景值、開邊界潮位值等。計(jì)算基面統(tǒng)一采用85國(guó)家高程基面。計(jì)算中，內(nèi)模時(shí)間步長(zhǎng)設(shè)為30 s，外模時(shí)間步長(zhǎng)設(shè)為6 s。垂向采用σ坐標(biāo)分為6層，垂直渦黏系數(shù)背景值取10-6(m2/s)，海底粗糙高度取值為0.001 m，最小底部拖曳力系數(shù)取0.001 5。模擬時(shí)間從2012年10月7日零時(shí)至2012年10月20日23時(shí)，海洋開邊界以潮位作為驅(qū)動(dòng)力，潮位值由東中國(guó)海潮波模型預(yù)報(bào)所得。

圖2 模型計(jì)算區(qū)域水深Fig.2 Bathymetry of study area

圖3 研究海域模型計(jì)算網(wǎng)格Fig.3 Model grid for study area

2.3 模型驗(yàn)證

為了確定模型精確性，引入了基于模型計(jì)算結(jié)果與實(shí)測(cè)資料一致性的預(yù)測(cè)能力系數(shù)，如下式：

(20)

該統(tǒng)計(jì)方法由Wilmott(1981年)[16]首次提出，Warner(2005年)[17]等和Li(2005年)[18]等都有使用。

圖 4 給出了油庫(kù)碼頭、開山島和翻身河口3個(gè)測(cè)站的潮位計(jì)算值與實(shí)測(cè)值的過(guò)程對(duì)比，模型計(jì)算的潮位值和相位結(jié)果均與實(shí)測(cè)值較為吻合，3個(gè)測(cè)站的潮位的預(yù)測(cè)能力系數(shù)都超過(guò)0.96(完全一致統(tǒng)計(jì)系數(shù)為1)。圖5和圖6分別給出了V1和V5測(cè)站大小潮表層和底層的流速和流向的計(jì)算值和實(shí)測(cè)值的過(guò)程對(duì)比。V1測(cè)站表、底層流速的預(yù)測(cè)能力系數(shù)分別為0.92和0.86，表、底層流向的預(yù)測(cè)能力系數(shù)為0.85和0.88。V5測(cè)站

圖4 計(jì)算潮位與實(shí)測(cè)潮位對(duì)比Fig.4 Comparisons of simulated and observed water surfer elevation

表、底層流速的預(yù)測(cè)能力系數(shù)都達(dá)到了0.96，表、底層流向的預(yù)測(cè)能力系數(shù)為0.94和0.89。可以看出，模型計(jì)算的流速計(jì)算值與實(shí)測(cè)值趨勢(shì)吻合，模擬結(jié)果的轉(zhuǎn)流時(shí)刻與實(shí)測(cè)值之間吻合較好,總的說(shuō)明該模型的模擬結(jié)果與實(shí)際的潮流整體上基本一致，符合研究海域的潮流分布特征。

2.4 潮流場(chǎng)特征分析

從大、小潮潮流矢量圖(圖7)可以看出：V1～V6垂線的潮流均主要表現(xiàn)為旋轉(zhuǎn)流。潮流從北往南逐漸增大，小潮時(shí)V1和V2的最大漲落潮流速在0.3 m/s左右，而V5的最大漲落潮流速分別達(dá)到了0.82和0.78 m/s。從流向上來(lái)看，研究海域中南部的V3、V4、V5、V6為典型的蘇北沿岸流，漲急時(shí)沿岸方向由南往北，落急時(shí)則改為由北往南。研究海域北部的V1、V2漲落急方向分別為向岸和離岸。大潮的流速地理位置特征與小潮類似，但流速明顯增大，V1的最大漲落潮流速分別達(dá)到了0.73和0.47 m/s，而V5的最大漲落潮流速分別達(dá)到了1.88和1.41 m/s，漲急流速明顯大于落急流速。

圖5 V1測(cè)站大小潮表層和底層流速、流向計(jì)算值與實(shí)測(cè)值對(duì)比Fig.5 Comparisons of simulated and observed velocity and direction at station V1

圖8(a)為大潮漲急流場(chǎng)分布圖，漲急時(shí)近岸以由南往北的沿岸流為主，灌河口以北海州灣區(qū)域漲急流向?yàn)槲髂舷?，流速大小分布?lái)看，海州灣內(nèi)流速較小，在0.4 m/s左右，由北往南流速逐漸增大，研究海域的流速在0.5 m/s左右，而射陽(yáng)河口海域流速達(dá)到了0.8 m/s。最大流速在射陽(yáng)河口以東外海，達(dá)到1.7 m/s左右。

圖8(b)為大潮落急流場(chǎng)分布圖，落潮時(shí)海州灣內(nèi)的流速方向?yàn)闁|北向，海州灣東北海域海水東-東北向流出計(jì)算區(qū)域，灌河口以北海域主要以東南向沿岸流為主，海州灣內(nèi)的流速大小在0.3 m/s左右，海州灣以南沿岸線走向，流速逐漸增大，最大流速在1.6 m/s左右。圖9小潮時(shí)漲落潮流速方向和大潮類似，漲急時(shí)海州灣內(nèi)的流速大小在0.3 m/s左右，落急時(shí)流速大小在0.2 m/s左右。最大漲落潮流速出現(xiàn)在射陽(yáng)河口以東外海，流速大小達(dá)到1.5 m/s左右。

圖7 大、小潮垂線平均流速矢量圖小潮大潮Fig.7 Depth averaged velocity neap tide spring tide

圖8 大潮表層流場(chǎng)Fig.8 Surface velocity field in spring tide

圖9 小潮表層流場(chǎng)Fig.9 Surface velocity field in neap tide

2.5 余流場(chǎng)特征

余流(ur,vr)可以由水平歐拉流速得到。由數(shù)值模型得到的瞬時(shí)歐拉流速(u,v)可以分解成周期項(xiàng)(up,vp) 和余流項(xiàng)(ur,vr):

(u,v)=(up,vp)+(ur,vr)

(21)

有2種計(jì)算余環(huán)流的方法。第一種方法是通過(guò)過(guò)濾或者時(shí)間平均瞬時(shí)流速，這種方法在很多文獻(xiàn)中都被采用，本文也將采用該方法。由于該方法易于使用，只需對(duì)模型計(jì)算結(jié)果進(jìn)行統(tǒng)計(jì)處理[7]。第二種方法需要通過(guò)求解余流演變方程，該方程通過(guò)對(duì)瞬時(shí)流速方程求取平均而得到。余流通過(guò)對(duì)瞬時(shí)流速求平均得到：

(22)

式中：T為潮周期，由于通過(guò)在潮周期內(nèi)求平均，大部分做超周期運(yùn)動(dòng)的變量通過(guò)求平均而消去。

從圖10和圖11來(lái)看，大潮時(shí)期表、底層近岸余流方向都是由北向南的沿岸流，灌河口以東海域，余流方向?yàn)闁|北向，海州灣以東外海余流方向東南轉(zhuǎn)東流向外海。南部海域的余流方向?yàn)槲鞅狈较颍c灌河口海域的海水混合后由東北向流出計(jì)算海域。海州灣內(nèi)的余流大小在0.01～0.03 m/s左右，海州灣以南的海域近岸余流大小在0.02～0.06 m/s左右，外海余流變大，最大達(dá)到0.15 m/s。底層余流明顯減小，海州灣內(nèi)的余流大小在0.01 m/s左右，近海區(qū)域余流在0.02～0.04 m/s左右，最大余流速達(dá)到了0.1 m/s。小潮時(shí)期，近岸余流方向還是以從北向南的沿岸流為主，海州灣外部海域靠近山東近岸流速方向?yàn)楸毕蛄鞒鲇?jì)算域，表層流速大小0.05～0.08 m/s，底層流速大小為0.03～0.05 m/s。計(jì)算區(qū)域北部和東南部邊界附近各有一個(gè)順時(shí)針漩渦，流速大小在0.1～0.15 m/s之間，由于潮流主要受潮波驅(qū)動(dòng)，余流漩渦形成的主要機(jī)理可能是潮波和地形及岸線的相互作用。南部海域海水從南部邊界進(jìn)入后分層北向和東北向兩股流，表層的北向流速大小在0.02～0.08 m/s，底層的北向流速大小在0.01～0.05 m/s；表層?xùn)|北向流速在0.06～0.15 m/s左右，底層?xùn)|北向流速在0.02～0.08 m/s左右。海州灣內(nèi)的表層余流流速大小在0.01～0.06 m/s，底層余流流速在0.01 m/s左右。

圖10 大潮余流流場(chǎng)Fig.10 Residual circulation in spring tide

圖11 小潮余流流場(chǎng)Fig.11 Residual circulation in neap tide

3 結(jié) 語(yǔ)

灌河口海域岸線地形復(fù)雜，該區(qū)域的水動(dòng)力過(guò)程具有較強(qiáng)的三維特性。本文以水平方向上采用三角形無(wú)結(jié)構(gòu)網(wǎng)格，垂直方向上采用σ坐標(biāo)FVCOM模式為基礎(chǔ)，建立了適應(yīng)于灌河口的三維潮流數(shù)值模型，能夠更好地?cái)M合了復(fù)雜岸形和海底地形。模型驗(yàn)證結(jié)果表明潮位、流速模擬值和實(shí)測(cè)值符合良好。

基于無(wú)結(jié)構(gòu)網(wǎng)格的FVCOM模型適用于蘇北灌河口海域，為下一步模擬灌河口泥沙和污染物變化過(guò)程提供了良好的水動(dòng)力場(chǎng)基礎(chǔ)。

□

參考文獻(xiàn)：

[1] 陳秀英，余乃旺，杭慶豐．平面二維水流數(shù)學(xué)模型在某項(xiàng)目防洪評(píng)價(jià)中的應(yīng)用[J]．人民珠江， 2011,(2):9-13．

[2] 黃家祥，殷 勇．灌河口潮灘重金屬累積特征及其對(duì)環(huán)境的意義[J]. 環(huán)境保護(hù)科學(xué)，2007，33(6):35-38．

[3] 董 佳，劉 勇，熊 偉. 灌河口水動(dòng)力條件對(duì)口門整治工程響應(yīng)分析[J].水運(yùn)工程, 2015,(10):125-131.

[4] Chen C S, Liu L, Beardsley R C. An unstructured grid, finite-volume, three-dimensional, primitive equation ocean model: application to coastal ocean and estuaries[J]. Journal of Atmospheric and Oceanic Technology, 2003,20:159-186.

[5] Isobe A, R C Beardsley. An estimate of the cross-frontal transport at the shelf break of the East China Sea with the Finite Volume Coastal Ocean Model[J]. Journal of Geophysical Research, 2006,111(C03012): 1-14.

[6] Weisberg R H, ZHENG Lianyuan. Circulation of Tampa Bay driven by buoyancy, tides, and winds, as simulated using a finite volume coastal ocean model[J]. Journal of Geophysical Research, 2006,111(C01005):1-20.

[7] ZHAO Liuzhi, CHEN Changsheng, Cowles G. Tidal flushing and eddy shedding in Mount Hope Bay and Narragansett Bay: an application of FVCOM[J]. Journal of Geophysical Research, 2006,111(C10015):1-16.

[8] CHEN Changsheng, HUANG Haosheng, Beardsley R C, et al. A finite-volume numerical approach for coastal ocean circulation studies: comparisons with finite difference models[J]. Journal of Geophysical Research, 2007,112(C03018):1-34.

[9] ZHENG Lianyuan. Hurricane storm surge simulations for Tampa Bay [J]. Estuaries and Coasts, 2006,29(6A):899-913.

[10] Lucy M Bricheno, Judith Wolf, Saiful Islam. Tidal intrusion within a mega delta: an unstructured grid modelling approach [J]. Estuarine, Coastal and Shelf Science, 2016,182:12-26.

[11] Michael Togneri, Matt Lewis, Simon Neill, et al. Comparison of ADCP observations and 3D model simulations of turbulence at a tidal energy site [J]. Renewable Energy, 2017,114:273-282.

[12] Mohammad Nabi Allahdadi, Chunyan Li. Effect of stratification on current hydrodynamics over Louisiana shelf during Hurricane Katrina [J]. Water Science and Engineering, 2017,10(2):154-165.

[13] 宋德海, 鮑獻(xiàn)文, 朱學(xué)明. 基于FVCOM的欽州灣三維潮流數(shù)值模擬 [J]. 熱帶海洋學(xué)報(bào), 2009,28(2):7-14.

[14] 李希彬, 鮑獻(xiàn)文, 馬 超, 等. 東海大堤對(duì)湛江灣水動(dòng)力環(huán)境影響的研究 [J]. 中國(guó)海洋大學(xué)學(xué)報(bào), 2009,39:287-296.

[15] Jiliang Xuan, Zhaoqing Yang, Daji Huang, et al. Tidal residual current and its role in the mean flow on the Changjiang Bank [J]. Journal of Marine Systems, 2016,154:66-81.

[16] Wilmott C J. On the validation of models[J]. Physical Geography, 1981,(2):184-194.

[17] Warner J C, Geyer W R, Kleczka J A. Numerical modeling of an estuary: a comprehensive skill assessment [J]. Journal of Geophysical Research, 2005,110 (C05001):1-13.

[18] LI Ming, ZHONG Liejun, Boicourt W C. Simulations of Chesapeake Bay estuary: sensitivity to turbulence mixing parameterizations and comparison with observations [J]. Journal of Geophysical Research, 2005,110(C12004):1-22.