蘇北灌河口海域三維水動力數(shù)值模擬

林偉波，魏愛泓，冒士鳳

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

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

1 三維有限體積模型及計算方法

1.1 控制方程

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

動量方程：

(1)

(2)

(3)

連續(xù)方程：

(4)

溫度和鹽度方程：

(5)

(6)

密度方程：

ρ=ρ(T,S)

(7)

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

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

(8)

(9)

1.2 邊界條件

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

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

(10)

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

(11)

阻力系數(shù)Cd由下面對數(shù)底邊界層計算值和常數(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)為表面凈熱通量，它等于短波輻射、長波輻射、感熱通量和潛熱通量之和；SW(x,y,0,t)為在海表面處短波輻射通量；cp為海水比熱系數(shù)。

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

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

(15)

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

(16)

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

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

(17)

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

(18)

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

側邊界處運動學、熱量和鹽度的邊界條件分別為：

(19)

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

1.3 數(shù)值方法

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

2 計算結果與分析

2.1 驗證資料說明

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

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

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

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

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

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

2.3 模型驗證

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

(20)

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

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

圖4 計算潮位與實測潮位對比Fig.4 Comparisons of simulated and observed water surfer elevation

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

2.4 潮流場特征分析

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

圖5 V1測站大小潮表層和底層流速、流向計算值與實測值對比Fig.5 Comparisons of simulated and observed velocity and direction at station V1

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

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

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

圖8 大潮表層流場Fig.8 Surface velocity field in spring tide

圖9 小潮表層流場Fig.9 Surface velocity field in neap tide

2.5 余流場特征

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

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

(21)

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

(22)

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

從圖10和圖11來看，大潮時期表、底層近岸余流方向都是由北向南的沿岸流，灌河口以東海域，余流方向為東北向，海州灣以東外海余流方向東南轉東流向外海。南部海域的余流方向為西北方向，與灌河口海域的海水混合后由東北向流出計算海域。海州灣內(nèi)的余流大小在0.01～0.03 m/s左右，海州灣以南的海域近岸余流大小在0.02～0.06 m/s左右，外海余流變大，最大達到0.15 m/s。底層余流明顯減小，海州灣內(nèi)的余流大小在0.01 m/s左右，近海區(qū)域余流在0.02～0.04 m/s左右，最大余流速達到了0.1 m/s。小潮時期，近岸余流方向還是以從北向南的沿岸流為主，海州灣外部海域靠近山東近岸流速方向為北向流出計算域，表層流速大小0.05～0.08 m/s，底層流速大小為0.03～0.05 m/s。計算區(qū)域北部和東南部邊界附近各有一個順時針漩渦，流速大小在0.1～0.15 m/s之間，由于潮流主要受潮波驅動，余流漩渦形成的主要機理可能是潮波和地形及岸線的相互作用。南部海域海水從南部邊界進入后分層北向和東北向兩股流，表層的北向流速大小在0.02～0.08 m/s，底層的北向流速大小在0.01～0.05 m/s；表層東北向流速在0.06～0.15 m/s左右，底層東北向流速在0.02～0.08 m/s左右。海州灣內(nèi)的表層余流流速大小在0.01～0.06 m/s，底層余流流速在0.01 m/s左右。

圖10 大潮余流流場Fig.10 Residual circulation in spring tide

圖11 小潮余流流場Fig.11 Residual circulation in neap tide

3 結 語

灌河口海域岸線地形復雜，該區(qū)域的水動力過程具有較強的三維特性。本文以水平方向上采用三角形無結構網(wǎng)格，垂直方向上采用σ坐標FVCOM模式為基礎，建立了適應于灌河口的三維潮流數(shù)值模型，能夠更好地擬合了復雜岸形和海底地形。模型驗證結果表明潮位、流速模擬值和實測值符合良好。

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

□

參考文獻：

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

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

[3] 董 佳，劉 勇，熊 偉. 灌河口水動力條件對口門整治工程響應分析[J].水運工程, 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] 宋德海, 鮑獻文, 朱學明. 基于FVCOM的欽州灣三維潮流數(shù)值模擬 [J]. 熱帶海洋學報, 2009,28(2):7-14.

[14] 李希彬, 鮑獻文, 馬 超, 等. 東海大堤對湛江灣水動力環(huán)境影響的研究 [J]. 中國海洋大學學報, 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.