郝子寧，張明亮，張云鵬，喬 洋

（大連海洋大學(xué) 海洋科技與環(huán)境學(xué)院，遼寧 大連116023）

大連到莊河海域位于黃海北部，其中包括老虎灘、大連灣和莊河部分海域。這部分海區(qū)為不規(guī)則的半日潮，其中M2分潮占主要成分，并且潮流呈往復(fù)流性質(zhì)。該海區(qū)的海洋資源條件優(yōu)越，海洋養(yǎng)殖生物種類繁多，資源量大。近年來，由于近海養(yǎng)殖業(yè)的迅速發(fā)展，莊河近海海域灘涂貝類、對蝦、人工流放增殖和浮筏養(yǎng)殖基地逐年增多，養(yǎng)殖水域的投餌、施肥、藥物使用等對近海域的水質(zhì)影響日益顯著，其近海養(yǎng)殖活動成了這片海域污染的主要來源之一。

隨著計算機(jī)和數(shù)值模擬技術(shù)的發(fā)展，采用海洋環(huán)境動力數(shù)學(xué)模型來研究海域的水動力和水環(huán)境狀況等問題受到越來越多的關(guān)注。如Moon基于ROMS（Regional Ocean Modeling System）的一個拉格朗日粒子追蹤實驗?zāi)M中國東北部黃海和日本海的巨型水母分布［1］。Mitaras應(yīng)用ROMS海洋模式，耦合拉格朗日粒子追蹤方法來研究海洋水體中物質(zhì)的輸運(yùn)結(jié)構(gòu)和路徑，并對污染物在海洋中的運(yùn)動狀態(tài)進(jìn)行預(yù)測［2］。Bilgili等基于二維模型耦合拉格朗日模塊研究粒子在美國新罕布什爾州河口的時空分布情況［3］。汪守東等采用水動力學(xué)模型POM（Princeton Ocean Model）模式模擬與預(yù)報了渤海海峽突發(fā)性溢油事故，采用Lagrangian追蹤法模擬油膜的輸移擴(kuò)散［4］。張越美等基于ECOM模型（Estuary Coastal Ocean Model）模擬渤海灣三維潮流場的時空分布特性［5］，王振華、龔煌及褚芹芹等［6-8］基于ECOMSED三維水動力模型，結(jié)合拉格朗日粒子追蹤技術(shù)，模擬長江口排污口、大連灣水域和膠州灣粒子的運(yùn)動軌跡及污染物輸運(yùn)過程。

近年來，結(jié)合有限差分的離散方法與有限元的網(wǎng)格易曲性等優(yōu)點(diǎn)的有限體積法被廣泛應(yīng)用，如Chen、Wei及Ge 等［9-11］采用三維海洋模型研究美國東南沿岸Satilla河、中國長江口和中國南海海水運(yùn)動特性，討論了網(wǎng)格的分辨率對其模擬精度的影響，王平等［12］基于三維潮流模型耦合拉格朗日粒子追蹤模塊模擬大連灣內(nèi)的水動力特性和特征粒子的運(yùn)動軌跡。這些模型是基于非結(jié)構(gòu)化網(wǎng)格的FVCOM海洋模式，可以更好地擬合不規(guī)則的岸線，在垂向坐標(biāo)計算上使用sigma坐標(biāo)變換，在垂向上使用相同的層數(shù)，使其在模擬不同水深、地形變化不劇烈時有很大優(yōu)勢，并且在模擬淺水區(qū)域也會得到較好的結(jié)果。目前，采用環(huán)境動力學(xué)模型來研究北黃海莊河養(yǎng)殖海域水動力和水體污染物運(yùn)動軌跡的較少，已有的研究是基于深度平均二維模型，這些模型對岸線的擬合精度及計算結(jié)果都不十分理想，并且該海域的水深變化顯著，因此應(yīng)更多地應(yīng)用三維模型來研究該海域養(yǎng)殖活動對其水質(zhì)的影響。

本文采用FVCOM海洋模型，耦合拉格朗日粒子追蹤方法，對北黃海莊河養(yǎng)殖海域進(jìn)行詳細(xì)的數(shù)值計算，基于拉格拉日拉格朗日粒子追蹤模塊模擬莊河養(yǎng)殖海域污染物的遷移規(guī)律，從潮流的角度分析該海域水體粒子的運(yùn)動過程。

1 數(shù)學(xué)模型

潮流模型為基于有限體積法的FVCOM數(shù)值模式，該模型結(jié)合了有限元體積法和有限差分法，便于離散計算海洋原始方程組，使其在理論和數(shù)值計算方面有很大優(yōu)勢。同時該模型兼顧擬合邊界和局部加密的優(yōu)點(diǎn)，潮汐的移動邊界用干濕判斷法處理。采用改進(jìn)的Mellor-Yamada 2.5 階紊流閉合模型及Smagorinsky湍流閉合法，使垂直和水平混合模型在物理和數(shù)學(xué)上閉合，不規(guī)則的底部邊界采用σ垂直變換來體現(xiàn)。

1.1 基本控制方程

控制方程采用σ坐標(biāo)系下三維連續(xù)、動量方程［13-14］：

式中 t為時間；u，v，ω為x，y，σ方向的速度分量；D為水深；η為基準(zhǔn)面到海表面距離。σ和z的關(guān)系式為σ=（z-η）/（H+η）；Km為垂直渦旋粘性系數(shù)；ρ0為參考密度，ρ為海水密度；g為重力加速度；f為科氏參數(shù)；DFx，DFy和D（Fq2，F(xiàn)q2l）為水平擴(kuò)散項，AM，Ah分別為水平和垂向紊動粘性系數(shù)；q2和l分別是湍動能和湍動的宏觀尺度；H為基準(zhǔn)面到海底的距離。

Km由修正的Mellor-Yamada（2.5）湍流閉合模型計算，計算公式為：

式中 Ps和Pb分別為剪切和浮力引起的湍流動能產(chǎn)生項；ε為湍流動能耗散率；Kq為湍流動能垂向渦粘擴(kuò)散系數(shù)；為近似壁面函數(shù)。

湍流動能混合長度方程由下列等式關(guān)系閉合：

Sm為穩(wěn)定函數(shù)。根據(jù)Mellor和Yamada，Sm由下列方程組確定：

粒子軌跡運(yùn)動方程可表達(dá)：

ω，w的關(guān)系：

1.2 邊界條件

在σ坐標(biāo)系下，海面σ=0，邊界條件為：

在海底σ=-1，邊界條件為：

式中τsx，τsy為海表面風(fēng)應(yīng)力；（u+v）為海底切應(yīng)力。

陸地邊界條件法向速度和通量為0，即vn=0。外海開邊界給定水位，用干濕邊界法處理開潮灘，注入海域的河流做源項處理。

2 莊河養(yǎng)殖海域的數(shù)值模擬

2.1 計算域與網(wǎng)格

海域水環(huán)境的演變過程主要受潮流的影響，也就是說，水體中的污染物是以潮流作為運(yùn)動的載體，在海水往復(fù)運(yùn)動的同時，其污染物會擴(kuò)散輸移至外海。莊河位于大連的東部，其海域有多個養(yǎng)殖區(qū)域，主要是底棲貝類生物。這類底棲貝類生物在養(yǎng)殖的過程中會產(chǎn)生、排泄一定的氮磷污染物，進(jìn)一步造成海域水體營養(yǎng)鹽富集，形成富營養(yǎng)化。鑒于莊河養(yǎng)殖海域獨(dú)特的特點(diǎn)，選擇該海域進(jìn)行污染物的擴(kuò)散輸移研究。該海域的計算區(qū)域范圍經(jīng)度為121.5°～124.1°，緯度計算范圍為38.4°～39.8°，包括老虎灘、大連灣、大窯灣和莊河部分海域。主要的匯入河流有5條，包括大洋河、沙河、碧流河、贊子河、登沙河，其入口給定流量邊界條件。

圖1為大連至莊河海域的平面網(wǎng)格示意圖，其中大窯灣、小長山島2個站點(diǎn)可以進(jìn)行潮位驗證，石城島附近的一個站點(diǎn)進(jìn)行潮流流速及流向驗證。海域地形來自中國人民解放軍海軍司令部航海保證部出版發(fā)行的北黃海海圖，利用SMS軟件提取地形，并生成無結(jié)構(gòu)化不規(guī)則三角形網(wǎng)格，水平方向采用可變分辨率策略，實現(xiàn)對近岸、海島及海底地形陡峭區(qū)域進(jìn)行局部網(wǎng)格加密，在近岸和海島附近，網(wǎng)格分辨率較高，例如石城島和小長山島養(yǎng)殖海域附近，在外海開邊界處及周邊的網(wǎng)格較為稀疏，如圖1所示。這樣劃分網(wǎng)格的目的在于精細(xì)網(wǎng)格能夠較好地復(fù)演重點(diǎn)區(qū)域的真實流態(tài)，外海粗網(wǎng)格可以節(jié)省計算時間，提高計算效率。整個計算海域包括13469個三角形網(wǎng)格，垂向設(shè)為10個σ層，外模時間步長設(shè)定為1s。該海域的潮位開邊界在計算域的下方，開邊界潮位過程線通過M2，S2，K1和O14個分潮計算。

圖1 計算域網(wǎng)格示意圖

2.2 監(jiān)測點(diǎn)的潮位對比

水動力模擬時間段選擇在2011年8月份，從8月16日至8月31日，包括一個完整的大潮和小潮過程。

圖2和圖3為大窯灣和小長山島潮位站模擬的潮位與實測數(shù)據(jù)的對比。由圖可見，無論是在大潮還是小潮時間段，大窯灣和長山島潮位的實測值和計算值吻合較好，誤差在0.2～0.5m范圍內(nèi)，總的來說模擬的誤差也在允許范圍內(nèi)。

圖2 大窯灣計算潮位和實測值比較（2011-08-16～2011-08-31）

圖3 小長山島計算潮位和實測值比較（2011-08-16～2011-08-31）

圖4～圖5為不同站點(diǎn)模擬的流速、流向與實測數(shù)據(jù)的對比，可見，模擬結(jié)果能夠反應(yīng)實際水流的運(yùn)動特征，潮流模擬的表層和底層的流向也基本一致，大潮時流速一般在0.3～0.6m/s之間，小潮時流速一般在0.1～0.3m/s之間。由于受海底地形影響，表層、底層流速差距為0.1m/s。

圖4 大潮期表層、底層的流速和流向計算值和實測值比較

圖5 小潮期表層、底層流速和流向計算值和實測值的比較

圖6為石城島海域周圍模擬的某一落潮時刻流場，最大流速出現(xiàn)在島嶼間的水道附近，在島嶼附近發(fā)生繞流，在島嶼或岸邊附近受地形影響流速明顯變小，整個區(qū)域潮流還具有明顯的往復(fù)特征。在黃海潮波的作用下，隨著潮位漲落，局部海域會出現(xiàn)周期性的偏渦區(qū)變化。

圖6 石城島周圍海域模擬的流場

2.3 模擬粒子運(yùn)動軌跡

模擬海中投放群粒子具有明顯的規(guī)律性，通過分析群粒子隨水體發(fā)生遷移和擴(kuò)散，可判斷污染物在該區(qū)域內(nèi)排放時對該水域整體水環(huán)境的影響。在石城島附近的養(yǎng)殖區(qū)域選取兩個投放點(diǎn)，通過拉格朗日法模擬群粒子的運(yùn)動軌跡來研究污染物的運(yùn)動規(guī)律。因模擬單個粒子的運(yùn)動規(guī)律存在一定的隨機(jī)性，為了更好地模擬污染物的運(yùn)動軌跡，采用群粒子為研究對象，粒子群運(yùn)移規(guī)律通過質(zhì)心軌跡來反映。

質(zhì)心坐標(biāo)的計算公式為：

式中 x（t），y（t），z（t）為粒子群質(zhì)心在某時刻的x，y，z軸的坐標(biāo)，xi（t），yi（t），zi（t）單個粒子在某時刻的x，y，z軸坐標(biāo)；n為模擬粒子總數(shù)。

粒子模擬時間段為2011年8月20日零時至2011年9月5日。2011年8月20日零時，向石城島東側(cè)和西側(cè)2個養(yǎng)殖海域（1#和2#）的水體中分別每層投放100個粒子，10層共放置1000個粒子，Begin為粒子投放點(diǎn)，End為粒子在結(jié)束時的位置。

粒子的運(yùn)動軌跡如圖7所示，通過模擬群粒子的輸移軌跡來研究污染物的運(yùn)動規(guī)律，發(fā)現(xiàn)粒子的運(yùn)動軌跡與釋放的位置有關(guān)。1#投放點(diǎn)位于石城島的東南側(cè)，2#投放點(diǎn)位于石城島的東北側(cè)。由圖7可以看出，通過15d追蹤，在1#投放點(diǎn)投放的粒子群在前幾天內(nèi)隨水流會做往復(fù)運(yùn)動，一段時間后粒子群質(zhì)心在潮流作用下向東南方向運(yùn)動，流向外海，粒子凈輸移距離較長，大約60km。該投放點(diǎn)的粒子群受潮流影響較大，易向外擴(kuò)散，不易發(fā)生滯留，能夠快速擴(kuò)散輸移至外海，計算結(jié)果表明該區(qū)的水體比較活躍，水交換能力較強(qiáng)。在2#投放點(diǎn)，在漲落潮流的作用下，粒子以螺旋形式沿著海岸線來回往復(fù)運(yùn)動，但總體向東側(cè)偏移，主要是落潮流驅(qū)動作用，粒子輸移的距離較短，大約20km，可以發(fā)現(xiàn)，該位置海域的水體交換較差，污染物擴(kuò)散、輸移能力較弱。同時也對石城島南側(cè)的部分海域進(jìn)行了粒子模擬，其粒子運(yùn)動朝向東南，一段時間后快速移向外海。

圖7 石城島周圍海域投放粒子的運(yùn)動軌跡

3 結(jié)語

北黃海海域附近島嶼眾多，并且海岸線極不規(guī)則，F(xiàn)VCOM模型基于有限體積法和非結(jié)構(gòu)化三角網(wǎng)格，其三角形網(wǎng)格能更好地擬合黃海北部海域的復(fù)雜不規(guī)則岸線和島嶼，同時有限體積法能保證每個網(wǎng)格單元的質(zhì)量守恒，進(jìn)而達(dá)到計算域整體質(zhì)量守恒。計算結(jié)果表明：

（1）污染物多伴隨落潮方向遷移，在靠近海岸方向遷移擴(kuò)散較慢，岸邊和島嶼附近的粒子滯留時間較長，所以在近岸水域大量的底棲貝類養(yǎng)殖將對岸邊的水環(huán)境帶來較大影響。

（2）石城島西南側(cè)海域的粒子凈輸移較長，該區(qū)的水體交換能力較好，而石城島東北側(cè)海域離岸較近，其粒子沿著岸線向東側(cè)運(yùn)動，其輸移距離較短，該區(qū)域的污染物擴(kuò)散輸移能力較弱，離岸線較遠(yuǎn)的石城島東側(cè)養(yǎng)殖海域污染物的輸移能力較強(qiáng)，遠(yuǎn)離海岸水域的貝類養(yǎng)殖水環(huán)境能夠在潮流的作用下進(jìn)行充分交換。

（3）各個投放點(diǎn)垂向上粒子運(yùn)動差別不大，主要原因是水深較淺，對環(huán)境影響較小。

［1］Moon J H，Pang I C，Yang J Y.Behavior of the giant jellyfish Nemopilema nomurai in the East China Sea and East/Japan Sea during the summer of 2005：A numerical model approach using a particle-tracking experiment［J］.Journal of Marine Systems，2010，80（1-2）：101-114.

［2］Mitarai S，Siegel D A，Watson J R，et al.Quantifying connectivity in the coastal ocean with application to the Southern California Bight［J］.Journal of Geophysical Research，2009，114，C10026：1-21.

［3］Bilgili A，Proehl J A，Lynch D R，et al.Estuary/ocean exchange and tidal mixing in Gulf of Maine Estuary：A Langrangian modeling study［J］.Estuarine，Coastal and Self science，2005，65：607-624.

［4］汪守東，沈永明，鄭永紅.海上溢油遷移轉(zhuǎn)化的雙層數(shù)學(xué)模型［J］.力學(xué)學(xué)報，2006（4）：452-461.

［5］張越美，孫英蘭.渤海灣三維變動邊界潮流數(shù)值模擬［J］.青島海洋大學(xué)學(xué)報（自然科學(xué)版），2002（3）：337-344.

［6］王振華，唐軍，李春良，等.長江口污染物運(yùn)動軌跡模擬［J］.海洋環(huán)境科學(xué)，2012（1）：1-5.

［7］龔煌，孫大鵬，孫志國.采用三維潮流模型分析人工島對海灣污染物運(yùn)動的影響 ［J］.水道港口，2013 （3）：263-271.

［8］褚芹芹，李磊，李培良.膠州灣潮流場的示蹤粒子數(shù)值模擬研究［J］.中國海洋大學(xué)學(xué)報（自然科學(xué)版），2010（11）：29-34.

［9］Chen C S，Liu H D，Beardsley R C.An unstructuredgrid，finite-volume，three-dimensional，primitive equations ocean model：application to coastal ocean and estuaries［J］.Journal of Atmospheric and Oceanic Technology，2003，20：159-186.

［10］Wei J，Zhang D，Xue P.et al.Coupling of a regional atmospheric model （RegCM3）and a regional oceanic model（FVCOM）over the maritime continent ［J］.Clim Dyn，2014，43：1575-1594.

［11］Ge J Z，Guo W Y，Ding P X.Hydrodynamic influence of proposed excavated-in harbor in the Hengsha Shoal of the Yangtze Estuary Ⅰ：Numerical Model and Validations［J］.Journal of East China Normal University （Natural Sc），2013，2013（4）：79-90.

［12］王平，張寧川.大連灣保守污染物遷移三維模型及應(yīng)用［J］.海洋通報，2013（3）：265-274.

［13］宋德海，鮑獻(xiàn)文，張少峰，等.基于FVCOM的廉州灣及周邊海域三維潮汐潮流數(shù)值模擬［J］.海洋通報，2012（2）：136-145.

［14］陳長勝.海洋生態(tài)系統(tǒng)動力學(xué)與模型［M］.北京：高等教育出版社，2003.