基于FVCOM的福建連江海域潮汐潮流數(shù)值模擬

于寒，劉桂梅，紀棋嚴，楊靜，朱學(xué)明，楊逸秋

（1.國家海洋環(huán)境預(yù)報中心，北京 100081；2.浙江海洋大學(xué)，浙江 舟山 316022；3.南方海洋科學(xué)與工程廣東省實驗室（珠海），廣東 珠海， 519000）

1 引言

福州市連江縣地處福建省東南沿海，位于閩江口北岸，是國務(wù)院首批沿海開發(fā)縣之一，縣內(nèi)有羅源灣、黃岐灣、定海灣3個灣。近年來連江縣大力發(fā)展海水養(yǎng)殖產(chǎn)業(yè)。連江海域位于福州市東南部海域，養(yǎng)殖密集且規(guī)模較大。連江海域也是福建省海洋災(zāi)害的多發(fā)海域，赤潮頻發(fā)、水環(huán)境溫度異常和缺氧等海洋生態(tài)災(zāi)害對養(yǎng)殖業(yè)的可持續(xù)發(fā)展造成影響。例如：2003年和2008年等多次發(fā)生養(yǎng)殖物因缺氧死亡的情況[1]；2012年5—6月在連江黃岐海域發(fā)生的米氏凱倫藻赤潮造成了巨大的經(jīng)濟損失[2]。連江海域?qū)儆诟＝ㄑ睾＃辈ㄏ到y(tǒng)由規(guī)則半日潮占主導(dǎo)[3]，對于以潮動力為主的羅源灣而言，潮汐和潮流是灣內(nèi)物質(zhì)運輸和水交換等水動力過程的物理基礎(chǔ)。利用數(shù)值模型研究該海域的基本水動力情況，掌握該海域潮汐潮流運動的主要特征，可為生態(tài)動力過程研究和污染物擴散分析提供基礎(chǔ)，對污染物防治和赤潮防控也具有重要的意義。

目前廣泛應(yīng)用的海洋模型包括：有限體積海岸海 洋 模 式（Finite-Volume Coastal Ocean Model，F(xiàn)VCOM）、區(qū)域海洋模式（Regional Ocean Modeling System，ROMS）、混合坐標大洋環(huán)流模式（Hybrid Coordinate Ocean Model，HYCOM）和普林斯頓海洋模型（Princeton Ocean Model，POM）等。FVCOM的非結(jié)構(gòu)三角形網(wǎng)格能較好地擬合近岸曲折復(fù)雜岸線，針對重點區(qū)域可進行局部加密，垂直方向采用σ坐標能較好地擬合海底地形。國內(nèi)外學(xué)者對近岸和海灣的潮汐潮流研究較多，例如：靖春生等[4]、岳云飛等[5]和林作梁等[6]基于FVCOM建立了不同海域較高分辨率的三維潮流數(shù)值模型，較好地再現(xiàn)了不同海域的潮流運動情況；CHEN等[7]模擬了緬因灣和新英格蘭大陸架的潮汐情況，分析了緬因灣半日潮和能量通量的近共振性質(zhì)以及新英格蘭大陸架上主導(dǎo)潮汐的復(fù)雜動力；CZIZEWESKI等[8]采用三維水動力模型對馬拉尼昂灣的潮汐環(huán)流進行數(shù)值模擬，研究了潮汐行為和相關(guān)能量的分布，發(fā)現(xiàn)河口地貌和淺層水深會大大增加各分潮的振幅和流速。一些學(xué)者也對羅源灣進行了研究，例如：杜伊[9]模擬了羅源灣不同深度處的潮流場，分析了灣內(nèi)外的水交換以及污染物輸運過程，結(jié)果表明納潮量的大小與水交換的強弱呈正相關(guān)關(guān)系；李真[10]利用FVCOM研究了羅源灣圍墾前后水動力的變化情況，并設(shè)計了圍墾前后以及不考慮灘涂情況下的納潮量和水交換率等水動力量，結(jié)果說明灘涂在水交換過程中具有重要作用。鑒于針對連江海域三維潮汐潮流特征的研究較少，且開展三維水質(zhì)模型需要以三維水動力模型為基礎(chǔ)，本文基于FVCOM建立了連江海域的三維潮汐潮流數(shù)值模型，模擬結(jié)果與觀測數(shù)據(jù)較為吻合，以此為基礎(chǔ)對研究區(qū)域的潮汐潮流特征進行分析，以期為該海域的生態(tài)動力過程以及污染物擴散研究提供水動力研究基礎(chǔ)。

2 模型建立

2.1 FVCOM模式介紹

FVCOM是陳長勝團隊開發(fā)的海洋模式[11]。模式主要特點為：①采用的有限體積法結(jié)合了有限元法易于擬合邊界、局部加密的特點以及有限差分法易于離散方程組的特點；②水平方向采用的三角形網(wǎng)格可更好地擬合岸線復(fù)雜曲折的海岸區(qū)域，垂直方向采用的σ坐標能夠較好地模擬海底地形；③采用干濕判別法能較好地模擬海水漲潮和落潮期間灘涂被淹沒或露出的情況[12]。連江海域海岸線曲折，地形復(fù)雜，F(xiàn)VCOM的非結(jié)構(gòu)三角形網(wǎng)格及σ坐標能滿足精細化模擬需求，還可針對重點關(guān)注區(qū)域進行加密。

2.2 模型配置

模型計算范圍為118.78°～121.34°E，24.33°～27.41°N。水平方向劃分為不重疊的非結(jié)構(gòu)三角形網(wǎng)格（見圖1a），在近岸和河道處加密處理，最小網(wǎng)格分辨率約為50 m，開邊界平均約為12 km，區(qū)域網(wǎng)格格點數(shù)為18 185個，單元數(shù)為33 124個；垂直方向采用σ坐標，分為等間隔的7層。模型采用內(nèi)外模態(tài)分離法求解，外模時間步長為1 s，內(nèi)模時間步長為5 s；水深數(shù)據(jù)應(yīng)用經(jīng)過國家海洋環(huán)境預(yù)報中心業(yè)務(wù)化海流數(shù)值模式測試的地形基礎(chǔ)資料，能覆蓋整個模型范圍，將地形數(shù)據(jù)插值到三角元節(jié)點上，形成計算網(wǎng)格節(jié)點的地形（見圖1b）。

圖1 計算區(qū)域模型配置Fig.1 Model configurations of calculation area

模型計算需要提供開邊界條件，開邊界條件的好壞對模擬結(jié)果有直接影響。兩種方法可以選擇開邊界，即提供開邊界節(jié)點的實時水位或調(diào)和常數(shù)。本文的開邊界采用OTPS（OSU Tidal Prediction Software）[13]的后報水位（中國海1/30°），包括了M2、S2、N2、K2、K1、O1、P1和Q18個主要分潮，可獲取與水動力有關(guān)的潮位邊界條件。

假設(shè)開始時海面靜止，即初始場所有三角元中心點的流速均為0，所有三角元節(jié)點的水位為0。采用正壓模式即只有潮位邊界驅(qū)動的模式。不考慮溫度和鹽度的時空變化，即假定兩者為常數(shù)，溫度為10℃，鹽度為32。不考慮河流的徑流輸入。時間為2013年1月1日00時—12月31日23時（北京時，下同）。

3 模擬結(jié)果驗證

為驗證模擬結(jié)果，潮汐觀測數(shù)據(jù)使用北礵（T1）、沙埕（T2）、三沙（T3）、秦嶼（T4）、長門（T5）和臺山（T6）6個臨時潮位觀測站以及東沃和崇武兩個驗潮站的水位數(shù)據(jù)。臨時潮位觀測站的數(shù)據(jù)分兩個時間段（見表1）。潮流觀測數(shù)據(jù)來自8個全潮觀測站點，站點位于寧德海域附近（見圖2），觀測時間和層次見表1。

圖2 潮位和潮流觀測站示意圖Fig.2 Schematic diagram of tidal level and sea current observation stations

表1 潮位、潮流觀測資料Tab.1 Observation data of tidal level and sea current

3.1 潮汐結(jié)果驗證

將模型的調(diào)和常數(shù)與驗潮站觀測數(shù)據(jù)的調(diào)和常數(shù)進行對比，結(jié)果見表2。整體而言，模擬得到的調(diào)和常數(shù)與觀測結(jié)果較為吻合。我們計算了M2、K1、S2和O1分潮振幅的平均絕對誤差分別為4.6 cm、2.5 cm、2.0 cm和1.1 cm，遲角的平均絕對誤差分別為3.7°、3.9°、6.5°和6.3°。從4個分潮振幅和遲角的相對誤差來看，個別站點振幅和遲角的相對誤差偏大，可能因為這些站點靠近岸邊，地形數(shù)據(jù)不準確，模型的岸線網(wǎng)格與真實岸線存在差距。

表2 M2、K1、S2和O1分潮調(diào)和常數(shù)觀測與模擬對比Tab.2 Observation and simulation comparison of harmonic constants of tidal components M2,K1,S2 and O1

對觀測水位進行調(diào)和分析可以得到各個觀測位置的調(diào)和常數(shù)，將觀測點的潮位與模式模擬結(jié)果進行對比（圖略）。8個驗潮站潮位觀測值與模式結(jié)果的平均絕對誤差分別為0.191 m、0.130 m、0.117 m、0.104 m、0.130 m、0.158 m、0.137 m和0.100 m。除個別觀測點在部分時段的潮位模擬值與觀測值存在較大差別外，整體上兩者的數(shù)值和相位變化都較為一致，模式基本可以模擬出潮位隨時間的變化情況。

3.2 潮流驗證結(jié)果

將潮流的模式結(jié)果與觀測數(shù)據(jù)進行對比。潮流觀測數(shù)據(jù)可分為大潮、中潮和小潮3個時段，垂向可分為6層（表層、0.2H、0.4H、0.6H、0.8H和底層）。在模式結(jié)果中，由于垂向上均勻分布的σ坐標可分為7層，流速位于σ層的中間層，因此可近似的將6個層次的流速與觀測數(shù)據(jù)進行對比。表3為8個站點在大潮、中潮和小潮期表層流速和流向觀測值與模擬值的平均絕對誤差。限于篇幅，本文只給出2#站點大潮期間表層、0.6H和底層的流速和流向?qū)Ρ惹闆r（見圖3）。從圖中可以看出，在部分觀測流速處于低值時，模擬結(jié)果偏大，模擬流速變化相位在部分時刻略早于觀測數(shù)據(jù)，除此之外，模擬值與觀測結(jié)果的變化趨勢基本一致。在大潮、中潮和小潮3個時段內(nèi)，8個站點流速模擬值與觀測值的平均絕對誤差的平均值分別為0.10 m/s、0.11 m/s和0.08 m/s，流向平均絕對誤差的平均值分別為24.27°，20.69°和23.47°，這表明模擬值可以較好地刻畫連江海域潮流的時空運動特征。

圖3 2#站點大潮期間表層、0.6H層和底層潮流的流速和流向?qū)Ρ菷ig.3 Comparison of sea current velocity and direction during spring tide at the surface,0.6H and bottom layers of Station 2#

表3 8個站點潮流表層流速、流向的平均絕對誤差Tab.3 Mean absolute errors of sea surface current velocity and direction of eight stations

4 潮汐、潮流和潮余流特征分析

4.1 潮汐特征分析

通過對計算區(qū)域水位的調(diào)和分析，我們獲得4個主要分潮的調(diào)和常數(shù)，以此計算潮汐類型的判別系數(shù)。具體公式及判別標準如下[14]：

式中：HK1、HO1、HM2分別為各個三角元節(jié)點上K1、O1和M2的振幅。計算得到海域所有節(jié)點潮汐類型的判別系數(shù)為0.11～0.40（見圖4a），這表明福建連江及其鄰近海域?qū)儆谝?guī)則半日潮類型，半日分潮占主導(dǎo)，即在1個太陰日內(nèi)出現(xiàn)兩個高潮和兩個低潮，潮汐的日不等現(xiàn)象不明顯。這種半日分潮占優(yōu)的現(xiàn)象在同潮圖中也有所體現(xiàn)（見圖5），主要表現(xiàn)為半日分潮的振幅普遍大于全日分潮。

對于規(guī)則半日潮海區(qū)，可以根據(jù)公式計算最大可能潮差。公式為[14]：

式中：HM2、HS2、HK1和HO1分別為各個三角元節(jié)點上M2、S2、K1和O1的振幅。通過計算可以得到每個網(wǎng)格點上的最大可能潮差（見圖4b）。由等潮差線可以看出，最大可能潮差與M2分潮等振幅線（見圖5）的分布一致，與岸線走向基本平行，潮差值由研究區(qū)域東北方向開邊界向近岸逐漸增加，從5.5 m逐漸增大到三沙灣和興化灣灣口的8.5 m，再向灣內(nèi)逐漸增加到9 m左右。由于淺水效應(yīng)和復(fù)雜岸線的影響，最大可能潮差越靠近岸邊越大，在各個海灣的頂部達到相對最大。

圖4 連江及其鄰近海域潮汐性質(zhì)Fig.4 Tidal properties in Lianjiang and its adjacent waters

圖5 連江及其鄰近海域主要分潮的同潮圖Fig.5 Co-tidal lines of main tidal components in Lianjiang and its adjacent waters

圖5為研究區(qū)域4個主要分潮的同潮圖，與福建沿岸4個主要分潮的同潮圖進行對比[14]，兩者的振幅、遲角的數(shù)值分布及形態(tài)都較為吻合。從等遲角線分布來看（見圖5），潮波主要沿著福建東南沿海自東北向西南方向傳播，與文獻一致[15]：潮波傳至東海陸架之后，多數(shù)直接到達浙江三門沿岸并形成高潮，后分為南北兩支，其中南支與從太平洋進入東海后向西南方向傳播的潮波匯合，一起傳到臺灣海峽。從等振幅線來看，受水深變淺和復(fù)雜岸線的影響，4個主要分潮的振幅都是越靠近岸邊數(shù)值越大，且等振幅線幾乎與岸線方向平行，4個分潮中M2分潮的比重最大，S2次之，K1和O1的比重相當(dāng)。

研究區(qū)域中4個主要分潮的振幅從開邊界到近岸逐漸增大，M2振幅增幅最大，超過80 cm，S2次之，增幅為30 cm，K1增幅為7 cm，O1增幅為5 cm。由此可見，相對于半日分潮的振幅變化，全日分潮的振幅增幅較小。由圖5可以看出，等振幅線和等遲角線在近岸海域基本垂直，表明研究海域內(nèi)的潮波主要為逆時針旋轉(zhuǎn)的駐波系統(tǒng)，此系統(tǒng)受臺灣海峽內(nèi)潮波系統(tǒng)的影響，在研究區(qū)域的近岸海灣（如羅源灣）中，等振幅線和等遲角線幾乎與岸線平行，表明海灣內(nèi)的潮波表現(xiàn)出前進波的特點。

4.2 潮流特征分析

本文采用類似于潮汐類型計算的方式計算潮流類型的判斷系數(shù)。具體公式及判別標準如下[14]：

式中：WK1、WO1和WM2分別為各個三角元上K1、O1和M2分潮流的表層最大流速。通過該公式計算得出（見圖6a），以平潭為界，平潭以北三角元上的潮流類型的判別系數(shù)小于0.5，屬于規(guī)則半日潮流，而平潭以南潮流類型判別系數(shù)多在0.5～2，屬于不規(guī)則半日潮流。這與莊樺[16]的研究結(jié)果相近，但本文的不規(guī)則半日潮流的分布更偏北一些。潮流的旋轉(zhuǎn)性質(zhì)由潮流橢圓要素中的旋轉(zhuǎn)率K決定，K是潮流橢圓的短半軸與長半軸的比值，即最小潮流流速與最大潮流流速的比值，當(dāng)K的絕對值大于0.25時，表現(xiàn)為旋轉(zhuǎn)式潮流，反之，表現(xiàn)為往復(fù)式潮流[17]。連江及其鄰近海域的K值等值線圖如圖6b所示，在靠近岸線以及海灣區(qū)域，K的絕對值多小于0.25，表現(xiàn)為往復(fù)流的形式；在興化灣、平潭、閩江以外海域以及東北海域內(nèi)，K值均大于0.25，表現(xiàn)為旋轉(zhuǎn)流的形式。

圖6 連江及其鄰近海域潮流性質(zhì)Fig.6 Tidal current properties in Lianjiang and its adjacent waters

對模式模擬的流速結(jié)果進行準調(diào)和分析，利用得到的潮流橢圓要素繪制4個主要分潮（M2、S2、K1、O1）的潮流橢圓圖（見圖7）。從圖中可以看出，灣內(nèi)海域潮流橢圓的短半軸幾乎等于0，長半軸方向幾乎與岸線平行，半日潮流的潮流橢圓長半軸大于全日潮流的長半軸。M2分潮潮流橢圓圖的長半軸最大為2 m/s，多為逆時針旋轉(zhuǎn)，僅在部分海灣以及狹窄水道里為順時針旋轉(zhuǎn)，且近岸以及海灣內(nèi)的橢圓長半軸方向大多平行于岸線；S2分潮的潮流橢圓圖的長半軸最大為0.6 m/s，與M2分潮類似，多為逆時針旋轉(zhuǎn)，且海灣內(nèi)長半軸多與岸線平行。從兩個半日分潮流的整體分布看，越靠近岸邊，潮流橢圓的長半軸越大，在進入灣口或狹窄水道后，長半軸達到最大。K1、O1分潮潮流橢圓的長半軸最大分別為0.19 m/s和0.17 m/s，多為逆時針旋轉(zhuǎn)，在興化灣口外海、平潭外側(cè)、福寧灣外海和三沙灣外海有一支順時針旋轉(zhuǎn)的潮流橢圓。

圖7 表層分潮的潮流橢圓圖Fig.7 Tidal current ellipses of surface tidal components

4.3 潮余流特征分析

在近岸海域，潮波在傳播過程中受到非線性作用產(chǎn)生了潮余流，所以水質(zhì)點經(jīng)過周期運動后，不會回到原來的位置。在計算過程中，將一個潮周期的潮流取平均可以得到歐拉余流，本文將多個周期的潮流進行平均，通過計算得到研究海域表層和底層的潮余流。羅源灣和三沙灣表層和底層的潮余流見圖8。一般而言，潮余流較小，在岸線劇烈彎曲的岬角處，當(dāng)潮流沿著海角流動時，岸線會因為曲率較大而出現(xiàn)較強的離心作用，使海水向外流去，產(chǎn)生離岸流。離岸流會使海角處的海面降低，兩側(cè)海水向海角附近補充，形成旋轉(zhuǎn)方向相反的兩個渦旋，離岸流右側(cè)為順時針方向旋轉(zhuǎn)渦環(huán)，左側(cè)為逆時針方向旋轉(zhuǎn)渦環(huán)[18]。

連江海域表層和底層潮余流的分布趨勢基本一致（見圖8）。底層潮余流明顯小于表層；潮余流流速整體較小，特別是在岸線平緩的海域以及水深較淺的開闊海域，而水深較大的狹窄水道處流速較大；北壁鄉(xiāng)為岬角地形，余流在海角處為離岸流，右側(cè)為順時針渦旋，左側(cè)為逆時針渦旋，流速明顯增大，可達20 cm/s；苔菉鎮(zhèn)也存在岬角地形，并出現(xiàn)明顯的離岸流以及右側(cè)順時針和左側(cè)逆時針的渦環(huán)結(jié)構(gòu)。

圖8 羅源灣和三沙灣海域潮余流分布Fig.8 Distribution of tidal residual current in Luoyuan Bay and Sansha Bay

5 總結(jié)

本文基于FVCOM三維海洋模型，水平方向采用非結(jié)構(gòu)三角形網(wǎng)格、垂直方向采用σ坐標，采用干濕判別技術(shù)以及有限體積方法，建立了福建連江海域的三維潮汐潮流數(shù)值模型。將潮位的模型調(diào)和常數(shù)與臨時驗潮站數(shù)據(jù)進行對比，結(jié)果較為吻合，潮位的平均絕對誤差小于20 cm；將潮流的模式結(jié)果與浮標數(shù)據(jù)進行對比，大潮、中潮、小潮流速的平均絕對誤差的平均值分別為0.10 m/s、0.11 m/s和0.08 m/s，流向平均絕對誤差的平均值分別為24.27°、20.69°和23.47°。對比結(jié)果表明模型能較好地反映連江及其附近海域潮汐潮流的時空分布特征。對模型結(jié)果進行分析，結(jié)論如下：

（1）連江海域的潮汐判別系數(shù)小于0.5，屬于規(guī)則半日潮類型；開闊海域的潮波主要表現(xiàn)為逆時針旋轉(zhuǎn)的駐波系統(tǒng)，羅源灣等海灣內(nèi)的潮波表現(xiàn)為前進波；最大可能潮差的走向與岸線基本平行，由研究區(qū)域東北方向開邊界向近岸逐漸增加，在各個海灣的頂部達到相對最大。

（2）連江海域的潮流類型以平潭為界，平潭以北潮流類型的判別系數(shù)小于0.5，屬于規(guī)則半日潮流，平潭以南的判別系數(shù)多在0.5～2之間，屬于不規(guī)則半日潮流；在靠近岸線以及海灣區(qū)域表現(xiàn)為往復(fù)流的形式，在興化灣、平潭、閩江以外海域以及東北海域表現(xiàn)為旋轉(zhuǎn)流的形式；潮流橢圓在海灣及水道內(nèi)的短軸幾乎為0，長軸方向與水道方向一致，4個主要分潮流的最大流速分別為2.0 m/s、0.6 m/s、0.19 m/s和0.17 m/s。

（3）表層、底層潮余流的分布形態(tài)基本一致，底層潮余流明顯小于表層，潮余流流速整體較小，在水深大的狹窄水道處流速較大，在苔菉鎮(zhèn)和北壁鄉(xiāng)的岬角地形處形成流速較大的離岸流，其左右兩側(cè)分別形成逆時針和順時針渦環(huán)。

本文考慮了正壓情況下的潮汐潮流特征分布，下一步將以此為水動力條件，模擬連江海域生態(tài)水質(zhì)等情況。