舟山綠色石化三維潮流特性數值研究

陸凡，沈良朵*，高郁，王晉寶

(1.浙江海洋大學 船舶與海運學院，浙江 舟山 316021)

1 引言

在經濟全球化快速發展的大背景下，世界各國對海洋資源的開發利用越來越重視。特別地，我國大部分沿海城市在港口碼頭、物流集散中心和綠色石化基地的建設投入逐步加大[1]。然而，部分以圍填海方案來獲取工業用地的建設方式由于缺少科學的指導依據，對周邊海域的生態環境、水動力條件和泥沙沖淤造成了不良影響[2]。基于此，目前通過采用數值模擬分析法以水工建筑物對周邊海洋環境的影響范圍為目標進行科學預測，能夠合理有效地解決如何綠色開發海洋資源的問題。

圍填海工程前能夠科學地模擬出建設海域相關的潮流運動形態差異和泥沙沖淤的顯著地區，對工程實施有著積極的推動作用。因此，一些國內外專家學者借助數值模擬進行了相關研究。任一晗等[3]基于FVCOM 三維水動力數值模型，選取1984 年、2010 年、2019 年3 個代表年份，探討圍墾工程影響下舟山群島海域潮流結構與潮能分布的時空變化狀況；Safavi等[4]采用MIKE3 平臺構建了烏魯米耶湖的三維數值模型，研究了受湖中堤道分隔影響湖中南北部水中鹽度的變化情況；竇明等[5]利用FVCOM 三維海洋數值模型，模擬嵊泗大洋山圍墾工程建設填海前后的流場及泥沙運動，對該工程附近海域潮流及各層泥沙分布變化進行了探討；劉嘉星和劉長根[6]利用驗證后的三維潮流數學模型ROMS 探討了天津港附近海域的潮流特征、潮流的分層特點和圍海造地引起的潮流場變化情況；劉金鵬等[7]以龍口市人工島為例，利用驗證合理的潮流場和波浪場數學模型，模擬分析了在極限風況條件下人工島群周圍海域水環境的變化特征及人工島岸線的波高特征值的變化規律。Zamani 和Koch[8]以伊朗西南部Maroon 水庫為研究對象，通過與實測數據進行比較的基礎上，對這兩種三維水動力模型的適用性進行了對比研究；Soltani等[9]通過Mike21，2011 Version 模型，模擬了海岸帶的波浪和水流，以確定奈班德海灣地區的泥沙輸移速率和方向；Elshemy 等[10]使用MIKE21 建模系統開發了曼扎拉湖的水動力模型和水質模型，研究了未來氣候變化對湖泊水動力和水質特征的影響。

雖然前人對圍墾工程的數值模擬研究已經取得不錯的結果，但大多數都是以FVCOM 三維模型或者MIKE 二維水動力模型對相關海域的水動力和溫鹽變化進行了細致研究，而采用MIKE3 和MIKE21 這兩種水動力模型對國內近海岸島嶼圍墾工程的對比性研究較少。基于此，本文以舟山綠色石化基地三期圍填海工程為研究對象，基于MIKE3 平臺構建了工程前后三維水動力數學模型，并考慮到MIKE21 中的二維模型也可以分析水平方向上的潮流運動情況，故對兩種模型進行了潮位驗證結果精度的比較分析。但結合到舟山群島海域水道縱橫、地形結構復雜的特性，水平方向上的二維研究無法呈現出三維模型展現出的流場空間結構變化，因此在三維模型下分析了采用“裁彎取直”和水道封堵方案下對潮流場特征的影響，并從表層、中間層和底層分別對工程前后的流速流向影響進行三維垂向剖析研究。

2 工程概況

本工程位于舟山群島和岱山以西的大、小魚山島海域，北瀕岱衢洋、南接灰鱉洋，周邊島嶼岸線曲折，港灣眾多。水下地形南北深、東西淺，海底高程基本在-65～-1 m，水下地形起伏明顯[11]。三期圍填海工程地處于121°55′11″～121°59′37″E 之間，其用海面積約為18.41 km2，成陸面積為17.09 km2。工程涉及海域潮流運動以往復流為主，而且基本為平行岸線方向的往復流。漲潮流主要為來自東海的前進波經由舟山群島之間的狹道進入杭州灣，其中岱山海域的漲潮流主流較為平順，整體上呈現由東南流向西北方向。落潮流則為杭州灣的落潮水流經舟山群島之間的狹道流入外海，落潮流基本與漲潮流一致，岱山周邊落潮流整體較為平順。海區內泥沙主要來源于長江口以及杭州灣外泄泥沙隨浙東沿岸流的輸移擴散，含沙量總體較大。

3 三維水動力模型

3.1 控制方程

三維模型使用的是分層網格，即在水平域中使用非結構網格，而在垂直域中使用結構化網格。垂直網格是分為sigma 坐標或組合的sigma/z 級坐標。本文采用的是考慮了自由曲面的sigma 坐標。使用sigma坐標最重要的優勢是它們能夠準確地表示水深測量，并在海床附近提供一定的分辨率[12]。采用單元中心的有限體積方法進行了原始方程的空間離散化，方程式如下：

式中，t為時間；x、y、z為笛卡爾坐標；η為表面高度；d為靜止水深；h=η+d為總水深；u、v、w為x、y、z方向的速度分量；f=2Ωsinφ，sinφ為科里力參數（Ω為角轉速，φ為地理緯度）；g為重力加速度；ρ為水密度；sxx、sxy、syx、syy為輻射應力張量的分量；νt為垂直湍流（或渦流）黏度；pa為大氣壓力；ρ0為水的參考密度；S為由于點源引起的排放量大小；us和νs為水被排放到環境水中的速度；Fu和Fv分別為x和y方向上的水平應力分項。

3.2 模型建立

3.2.1 計算區域及網格

考慮到工程區域附近漲落潮的流路經杭州灣和部分東海海域，因此所選取區域范圍遠大于三期工程計算區域，北界至32°N 左右，南界至28.5°N，東到125°E 附近，潮流界頂點江陰作為長江上游邊界，澉浦潮位站附近作為杭州灣上游邊界，計算區域位于120°～126°E 之間。三維模型采用sigma 網格，垂向均勻分10 層，工程附近海域網格進行加密，網格最小尺度為18 m，工程外海區域，網格距為8 000 m，不同尺度網格之間通過設置實現平滑過渡，共計52 007 個節點，99 195 個單元。對于二維模型，共計41 041 個節點，78 473 個單元。模型計算區域及水深如圖1 所示。

圖1 模型區域和水深Fig.1 Computational domain and water depth

3.2.2 模型參數設置

模型計算時間步長設置根據克朗數CFL 條件進行動態調整，在確保模型計算穩定進行下，平均時間步長為2 s，最短時間步長為0.01 s。干濕水深判別設置為干水深0.005 m，淹沒水深0.05 m，濕水深0.1 m。這有利于保護模型穩定運行。水平渦黏系數經率定取0.28，垂直渦黏系數使用對數公式確定。通過粗糙高度來表征底床糙率，率定后取0.003 m。邊界采用水位控制，即用潮位預報的方法得到開邊界條件。外海開邊界潮位根據由16 個主要分潮（M2、S2、N2、K2、K1、O1、P1、Q1、MU2、NU2、T2、L2、2N2、J1、M1和OO1）的調和常數，以開邊界潮位預報公式計算得到：

式中，A0為平均海面；Fi和 (ν0+u)i為天文要素；Hi和 gi為調和常數。在閉邊界取流速的法向導數為0，在潮灘區采用漫灘邊界處理。

4 模型驗證

由于本工程周邊缺少實測數據，故本次模型驗證采用了距魚山島較近的小洋山海域，由中交第三航務工程勘察設計院有限公司在水文測驗期間的潮流實測數據。其中，分大中小潮對臨時測站T1 和T6 進行潮位驗證（大潮測流時間為2015 年5 月4 日至5 日；中潮測流時間為2015 年5 月7 日至8 日；小潮測流時間為2015 年5 月12 日至13 日）；同樣地，對定點測站N1、N14 以及走航測站WP05 流速流向數據進行潮流分層驗證，各測站位置見圖2。

圖2 驗證點位置Fig.2 Verification point location

采用MIKE3 的水動力模型計算數據與實測數據進行驗證，潮位驗證結果如圖3 所示。圖中顯示，從整體上看T1 測站和T6 測站潮位擬合較好，其中小潮階段擬合效果最好，相對誤差結果在±10%左右，個別高潮位處相對誤差稍大。總體來看，三維模型模擬的潮位驗證是符合要求的。

圖3 潮位驗證Fig.3 Tide level verification

對潮流測站N1、N14 和WP05 分為表層（第10層）、中層（第5 層）和底層（第1 層）進行流速流向驗證。在分層驗證下，可以清楚看到表層小潮階段的流速擬合結果是相對較大的，其他各層擬合結果都是極好的。表層模擬結果如圖4 所示，各測站僅在高潮位由于流量加大的影響導致流速擬合結果有0.5 m/s左右的誤差，其他時間段擬合相對誤差控制在8%左右，表層流向模擬結果令人滿意。中間層模擬結果如圖5 所示，各測站除高潮位處流速存在0.15～0.3 m/s 的誤差，其余各時間點相對誤差在10%以內，流向擬合結果較滿意。底層模擬結果如圖6 所示，大潮和中潮階段的模擬流速僅在N1 測站的高潮位時有0.35 m/s 的誤差，其余擬合較好，小潮階段除在個別時間點有0.45 m/s 的誤差外，其余各點存在0.15～0.25 m/s 的誤差范圍；底層流向模擬結果較好，流向相對誤差在7%以內，極個別點處由于測量過程存在失誤而出現的誤差，不予考慮。

圖4 表層流速流向驗證Fig.4 Verification of flow velocity and direction in the surface layer

圖5 中間層流速流向驗證Fig.5 Verification of flow velocity and direction in the middle layer

圖6 底層流速流向驗證Fig.6 Verification of flow velocity and direction in the bottom layer

為對模型驗證結果有更加科學性的認識，此處通過Willmott[13]所提出的skill 方程進行評估，表達式如下：

式中，Pi是模擬值；Oi是 實測值；是實測平均值。通過d值判斷模型效率，d＞0.65 時，模型效率極好；0.5＜d≤ 0.65，表示模型效率非常好；0.2＜d≤ 0.5，表示模型效率好；d≤ 0.2，表示模型效率差。

結果表明，N1 和WP05 測站各層的大中小潮期間流速d值在0.80～0.97 之間，流向d值區間為0.73～0.99，其模型效率都是極好的。限于篇幅，僅具體列出了N14 測點各層的模型效率值，具體見表1。

表1 模型效率系數Table 1 The index of agreement

5 三維結果分析

5.1 MIKE3 與MIKE21 對比分析

對MIKE3 與MIKE21 兩種擬合效果進行對比分析，分別將這兩種水動力模型的計算數據與實測數據進行驗證，潮位對比驗證結果如圖7 所示，由圖可知：兩種模型的潮位計算值與實測值的吻合度較高，特別是在小潮期間，其最大誤差僅有0.15 m 左右。大潮和中潮期間，三維數值結果與實測值在個別高潮處有0.5 m 左右的誤差，總體上看，在將二維模型中的初始條件同樣地設置于三維模型的情況下，二維數值結果更接近于實測值。

圖7 潮位對比Fig.7 Tide level comparison

考慮到理論方面上MIKE3 中間層的潮流模擬結果與MIKE21 的平均流速流向的結果相似，故在此處進行了流速流向的對比驗證。中間層模擬結果對比如圖8 所示，圖中各測站流向驗證結果良好。本文選擇的水深0.6H的實測數據與垂向平均流速基本保持一致，在大潮和中潮階段兩種模型擬合結果較好，最大相對誤差在20%左右，相較而言，小潮階段的高潮位處誤差稍大，相對誤差在25%左右，其余各時間點兩者模擬結果均良好。通過這4 組數據的對比發現：二維模型的大部分計算值更加接近于實測值，尤其是在漲急時刻，原因可能是二維模型是取值于垂向水深平均的流速進行模擬，而本文三維模型未考慮風場應力的影響，表層流速計算誤差會對三維流速計算結果有影響。總體上看，三維模型中間層的模擬結果與二維計算結果是最接近的，且中間層擬合效果最為出色，模擬精度基本滿足《海岸與河口潮流泥沙數值模擬技規程》的要求，模型可用于工程后的預測等各項工作[14-15]。

5.2 工程前后潮流場分析

從計算的潮流場分析得出，漲急時刻外海潮波沿東南向西北方向前進。在大魚山的東南側受島嶼阻擋，潮流分成兩股，分別沿著大魚山島東南側向南北流入杭州灣。大魚山島區域漲潮時刻受峽道效應，南岸流速大于北岸，落潮時刻則北岸較大，工程南北兩端的磯頭岸線對漲落潮流起到了挑流的作用，而且兩側都為深槽，因此漲落潮流速相對較大，但整體還較為平整。其中，漲潮時平均流向為230°～300°，落潮時平均流向為70°～110°，漲落潮平均流向約為180°，基本符合往復流性質。

考慮到圍墾工程所處的大小魚山島周邊范圍內島嶼錯落，水道縱橫，潮流運動情況復雜，故對工程前后計算域內大潮漲落時刻的表層流場變化進行具體的對比分析，同時給出了工程后流速的等值線圖以更直觀看出流速大小變化，如圖9 所示。整體上看，此圍填海工程前后對工程周邊海域的潮流運動沒有較大影響，是明顯的往復流形式。在大潮漲急階段，由于受圍墾工程“截彎取直”方法的影響，處于工程東北側的岸線外伸，東北側岸線存在一個向內的45°轉角，受岸線角度和洋流反射作用，較于工程前平行于岸線的潮流運動方向，此處產生一個直徑約為700 m 的逆時針旋渦；魚山島的西北側和東南側受圍填影響，潮流運動遇到外擴海岸線，此處潮流運動方向變化呈現出30°轉角；大范圍流場前后變化不大。大潮落急階段，工程西北角處的一處小島嶼圍填時與主島進行了連接處理，將此處原先經狹道效應的高速流與經岸線折射過來的洋流形成的旋轉流影響削弱，故旋轉流轉變成了與岸線彎口一致的運動形式；魚山島東北側、南側以及西側落潮流受外延海岸線的挑流作用，工程后流場運動均出現30°～45°的轉角變化，這表明處于迎流面的岸線變化對流場的改變起主導作用。

圖9 表層潮流場對比（上：漲急 下：落急）Fig.9 The comparison of tidal current field on the surface layer (up:rising;down:falling)

5.3 垂直流場分析

為更立體的表現出工程區域的流速流向情況，本節同時選取了N1、N14 和WP05 這3 個典型的測站數據進行垂直結構的分析。由圖10 分析得：N1 測站的流速在垂直方向上的實測值和漲落急階段的計算值各自都較接近，并在中間層附近的水深0.6H處僅有4.6%左右的相對誤差，計算與實測的最大誤差出現在表層附近，相對誤差在-25%左右，漲落急各分層流速較為一致，整體上都呈現出流速隨水深增加而減小，且幅度較大，可能是受到底部海床的摩阻力較大的影響；對于各分層流向趨于一致，且與實測數據擬合較好。N14 測站處漲落急時刻的流速計算值與實測值最大誤差出現在表層，其相對誤差為-28.01%，同時此測站的漲急流速計算值大于落急時刻計算值，這可能是由于此處網格布置較疏的影響，加上缺少風場實測資料，未考慮風場的原因所導致。WP05 測站可能處在洋山港與殼子山之間的狹窄水道里，由寬闊海域進入該水道，過水斷面突然減小，故同樣是在表層處的流速與實測數據有-30%的相對誤差，而0.6H以下水深的計算流速與實測值的誤差范圍為13～22 cm/s，且漲急時刻流速小于落急流速，同樣地各分層的流向是近乎一致的。通過對垂直流場的分析，可以清楚地得到流速與水深呈負相關，流向受水深影響近乎不計。而針對表層擬合較大問題，初步考慮是受到了風場和波浪對潮流的驅動作用。王世澎等[16]在二維情況下波浪對潮流場的數值分析中指出，對于流速過程，波浪在沒有破碎的情況下對流速過程的影響比對潮位過程的影響大。考慮到圍填海工程后，周圍海域受局部挑流作用，加之波浪輻射應力影響，波浪會在近岸破碎產生沿岸流疊加到原有潮流場上，對表層流速影響較大。針對此問題，將會在以后的研究中通過波浪耦合進行率定工作。

圖10 測站垂向流速流向分布Fig.10 Flow velocity and direction in vertical direction

6 小結

本文以舟山綠色石化三期工程為背景，采用MIKE3建立寧波舟山及其周邊附近海域三維水動力模型。通過與實測數據進行潮流驗證，對比了二維模型與三維模型驗證結果的區別，并主要對圍墾工程前后的水動力變化情況進行了詳細分析。最終得到如下幾點結論：

（1）以MIKE21 FM 的模型參數設置于MIKE3 中，二維潮位驗證效果更好一點，但兩種模型的驗證都是符合要求的；兩種模型流速擬合的最小相對誤差都控制在10%左右，通過4 組數據的對比，更直觀地驗證了二維數模數據是取值于垂向平均流速的理論結果，而三維中間層的模擬結果與二維計算結果則是非常接近。

（2）圍墾工程前后并未大面積改變大魚山島附近的潮流運動，僅在工程東北側、東南側和西北角由于圍墾工程將小島嶼之間的水道連接起來，并在“截彎取直”的圍墾方案影響下，使得局部區域潮流流態發生改變。

（3）以MIKE3 對工程附近垂向分布上的流速流向進行了分層結構分析，整體上流速都呈現出隨水深的增加而減小，可能由于未考慮風場和波浪場的影響，導致部分測站表層計算數據與實測數據誤差稍大，對于流向各層十分一致。