珠江河口潮流的數值計算研究

王光宇，宏波

(華南理工大學土木與交通學院，廣東廣州510641)

珠江河口水系復雜，網河區水網稠密，河道縱橫交錯，水流相互貫通，具有典型的潮汐網河水特征，動力過程變化多端，存在年、季、月、周的周期變化。在徑流、潮汐、季風、沿岸流和南海暖流等因素的綜合作用下，珠江河口存在復雜多變的水動力條件[1]。分析珠江河口的潮流作用，有助于進一步認識珠江河口的水動力過程。丁芮等[2]采用無結構網格三維有限體積海洋模式FVCOM，建立了珠江口及鄰近海域的三維正壓高分辨數值模型，發現珠江口海域潮汐為不正規半日潮，M2分潮占主導地位，分潮在淺海陸架區以旋轉流為主。Mao等人[3]通過對實測資料計算分析發現珠江口水域的潮差在1～1.7 m之間，M2分潮為主要的潮流成分。劉歡等[4]利用ECOMSED-3D數學模型，計算了珠江河口的能量傳播和能量耗散特征。鄭舒等[5]利用ECOMSED模型研究發現珠江口以M2分潮為主，潮振幅由外海向口門傳播過程中逐漸增大。包蕓[6]采用Backhaus三維斜壓模型模擬西南風場對珠江口近岸海域的影響，發現風可以明顯改變水流運動的大小。許煒銘等[7]對珠江河口進行了枯季二維水動力整體模擬計算，發現漲落潮流在各口門不同河道之間有明顯的非同步性。此外，陳翔等[8]構建了贛江的SWAT分布式水文模型，金浩宇等[9]利用SWAT模型模擬了尼洋河流域徑流過程，姜容等[10]建立了西江河口河網水質數學模型，模型驗證結果良好，為河網地區的相關研究提供了一定的幫助。

潮流是潮汐河口地區主要水動力條件之一，是河口地區最基本的物質運動。河口地區的島礁建設、橋梁建設等海上工程都會影響當地的潮流流動，并進一步影響泥沙、鹽分、污染物等物質的運輸。建立河口區域的數值模型，可以很好地預估海上工程對周遭流場的影響，有助于優化工程的設計方案及評估工程結構物對海洋環境的影響。趙強[11]等計算了圍填海工程對南黃海輻射沙脊群海域的潮流和含沙量變化。張楊[12]模擬計算了人工島對湯河河口動力過程的影響，發現工程前后水流流態發生明顯變化。季榮耀等[13]研究了港珠澳大橋人工島對水動力的影響，發現人工島兩端繞流區流速增大14%，背水面回流區最大減幅約27%，并存在強紊動小尺度回流。李雨[14]利用Mike軟件驗證了斯里蘭卡漢班托塔港人工島的設計合理性。王李吉等[15]研究發現瓊州海峽人工島建設會對工程區域產生一定的影響。本研究以EFDC模型為基礎，建立了珠江口水域的潮流數值模型，經模型驗證，此模型的計算結果與實測資料吻合良好，為河口地區的潮流數值研究提供了一種模型工具。

1 模型介紹

EFDC水動力學方程假定垂向靜壓，水平方向采用曲線正交坐標系，水質方向采用δ坐標變換，沿重力方向分層，求解三維紊動黏性方程，水平邊界擬合正交曲線坐標系和垂向δ坐標系下控制方程。

動量方程：

?t(mHu)+?x(myHuu)+?y(mxHvu)+

?z(mwu)-(mf+v?xmy-u?ymx)Hv=

-myH?x(gξ+p)-my(?xh-z?xH) ?zp+

?z(mH-1Av?zu)+Qu

(1)

?t(mHu)+?x(myHuv)+?y(mxHvv)+

?z(mwv)-(mf+v?xmy-u?ymx)Hu=

-mxH?y(gξ+p)-mx(?yh-z?yH)?zp+

(mH-1Av?zv)+Qv

(2)

?t(mξ)=-gH(ρ-ρ0)ρ0-1

(3)

連續方程：

?t(mξ)+?x(myHu)+?y(mxHv)+?z(mw)=0

(4)

(5)

狀態方程：

ρ=ρ(P,Sa,T)

(6)

溫度輸運方程：

?t(mHT)+?x(myHuT)+?y(mxHvT)+

?z(mwT)=?z(mH-1Kv?zT)+QT

(7)

鹽度輸運方程：

?t(mHSa)+?x(myHuSa)+?y(mxHvSa)+?z(mwSa)=?z(mH-1Kv?zSa)+QT

(8)

式中u、v、w——邊界擬合正交曲線坐標x、y、z方向上的速度分量；t——時間；mx、my——水平坐標的變換因子；m——度量張量對角元素的平方根，m=mx·my；Av——垂向紊動黏滯系數；Kv——垂向紊動擴散系數；f——科里奧利系數；ρ——混合密度；ρ0——參考密度；H——總水深；h——未擾動的z坐標原點以下的水深；p——壓力；Sa——鹽度；T——溫度；ξ——自由地勢能；Qu、Qv——動量在x和y方向的源匯項；QT——溫度的源匯項。

2 模型建立

珠江河口水流的重要特點之一，是其縱橫交錯的河網之間、河網與八大口門之間相互貫通、交叉影響；同時，在上游徑流和潮汐動力的共同作用下上下游往復流動，在各河道和河口內變化非常復雜。模型的計算區域包括河網區、河口灣區和近海水域，計算范圍為21°30′～22°49′N，112°30′～115°30′E。網格數為263×246，x方向網格最小長度68.89 m，最大長度4 846 m；y方向網格最小長度39.36 m，最大長度5 845 m(圖1)。上游邊界由實測水位或流量給出，外海邊界由調和常數推出。計算采用三維正壓模式，垂向分10層。潮汐諧波數據由俄勒岡州立大學潮汐預測軟件(OTPS)算得，并用于模型開邊界的驅動。開邊界的溫度和鹽度數據取自WOA13[16]，風力數據來源于香港機場觀測站采集的數據。河流徑流量數據采用珠江3條主要支流的上游各站點，即石角站(北江)、高要站(西江)和博羅站(東江)的觀測數據。

圖1 模型計算區域及網格劃分

3 模型驗證

取模型運算1 a(360 d)的計算結果與燈籠山、黃埔、九州港、內伶仃、桂山島、三灶、MO1共7個驗潮站1 a的實測水位數據進行對比驗證，模型時間步長為5 160，底摩擦系數按公式，其中卡門參數取0.4，糙率高度為0.003，為模型底層到海底的高度。站位分布情況見圖2。

為方便顯示，取燈籠山、黃埔、九州港、內伶仃、桂山島驗潮站前30 d的驗證結果進行表示。取三灶、MO1驗潮站后30 d的驗證結果進行表示。模型的具體驗證結果見圖3、4，模型的計算結果與實測值吻合良好，除黃埔和內伶仃驗潮站的水位誤差略大為0.2 m，其他各站均方根誤差在0.15 m左右。

流速驗證采用2007年MO1站的全年觀測資料，并分別對洪季(7月，模型計算第181—210天)和枯季(12月，模型驗證第331—360天)，水深為4 m和5 m的中層水體的流速計算結果進行了驗證，驗證結果見圖5、6。

“●”表示驗潮站站點，“△”表示后文潮流分析的站點圖2 站點分布位置

a) 燈籠山

b) 黃埔

c) 九州港圖3 燈籠山、黃埔、九州港、內伶仃、桂山島站水位驗證結果

d) 內伶仃

e) 桂山島續圖3 燈籠山、黃埔、九州港、內伶仃、桂山島站水位驗證結果

a) 三灶

b) MO1圖4 三灶、MO1站水位驗證結果

a) 水深為4 m的流速

b) 水深為5 m的流速

圖5MO1站洪季流速矢量驗證結果

a) 水深為4 m的流速

b) 水深為5 m的流速

圖6MO1站枯季流速矢量驗證結果

4 潮流分析

為了描述該區域的潮流分布特征，分別選取了7個具有代表性的站點，其中1、2號站點位于伶仃洋海區，3號站點位于內伶仃島以南，4號站點位于珠海橫琴島以南，5號站點位于萬山群島以西，6號站點位于三門島側灣海區，7號站點位于外伶仃島至擔桿列島之間海區。站點分布位置m見圖2中“△”。丁芮等[2]指出，珠江河口水域的主要潮流成分為M2、K1、S2、O1分潮。經過計算發現O1和S2分潮的影響相對較弱，為精簡表述，故取M2分潮和K1分潮作為主要的分析對象。

4.1 潮流性質

潮流性質主要根據潮型數Z=(WO1+WK1)/WM2來判別。若Z≤0.5則為正規半日潮，若0.54.0則為正規全日潮。除5號站位表層外，珠江河口水域潮流的潮型數為1.3～1.8(表1)，屬不正規半日潮流。這與前人發現的珠江河口水域主要以不規則半日潮為主的結論相一致[2-3,5]。

表1 各站點表層與底層潮型數

4.2 潮流類型

根據橢圓率(橢圓短軸與長軸之比)可以判斷潮流的類型。當橢圓率為0時，表現為嚴格的往復流，當橢圓率為1時，表現為理想的旋轉流。除淺海陸架海域6、7號站點橢圓率均大于0.5，1、5號站點的橢圓率也接近0.5，而其余站點的橢圓率均小于0.3，表明珠江河口水域主要以往復流為主，這主要是因為徑流下泄和狹窄的河口地形等因素的影響，使得珠江河口內部不能形成回轉流，只能形成在2個相反方向上流動的往復流，而寬闊的淺海陸架海域受地形的限制較小，在多個潮波的干涉作用和地轉偏向力的影響下，潮流類型以旋轉流為主(表2)。

表2 各站點M2和K1分潮表層與底層潮型數

4.3 潮流流速的季節變化

分別統計珠江河口水域底層和表層冬(1—2月)夏(7—8月)兩季的平均潮流流速。從底層的潮流平均流速季節變化(圖7a)中可以看出，除5號站點外，模型區域冬夏兩季潮流的平均流速的大小基本保持一致，不發生明顯變化，流場性質較為穩定，潮流方向均以北向為主。冬夏兩季最大的平均流速均發生在2號站點。對比與近底層的潮流流速，模型區域表層的流速存在一定的季節變化特性，這與徑流和季風的季節變化有關。從夏季到冬季，1、6號站流速明顯變小，3、4、5號站的流速變大，2、7號站點流速大小較為穩定。模型區域的流向發生明顯的西向偏轉，由夏季的偏東方向變為冬季的偏西方向。其中6、7號站點的流速方向發生180°偏轉，由東南向轉為西北向(圖7b)。相較于穩定的底層潮流流速，表層潮流流速的季節性變化應當與珠江口徑流下泄、季風風向以及外海潮波的季節性變動有關。

a) 底層

b) 表層圖7 潮流流速的季節變化

4.4 潮流橢圓的季節變化

選取夏季與冬季的潮流數據進行對比，并分別繪制M2分潮和K1分潮在夏季與冬季的潮流橢圓。在圖8中以紅色線表示夏季潮流橢圓，以藍色線為冬季的潮流橢圓。為了方便表示，對表層潮流橢圓縮小20倍表示，對底層潮流橢圓縮小50倍表示。

4.4.1M2分潮

M2分潮近底層潮流類型不隨季節發生變化，1號站點近底潮流表現為帶有一定旋轉性質的往復流，2—4號站點均表現為往復流；5—7號站點表現為旋轉流。各站夏冬兩季潮流橢圓的大小基本不變,但旋向相反，1—5號站的潮流橢圓旋轉方向由夏季的順時針變為冬季的逆時針，6、7號站則由夏季的逆時針旋轉變為冬季的順時針旋轉(圖8)。在潮流流速分布方面，伶仃洋海區的2、3號站點、狹長珠海海區的4號站點和近陸架寬闊海區的5—7號冬季與夏季潮流流速十分接近，僅大濠島以南海區的1號站點是冬季潮流流速大于夏季的，冬季最大近底流速出現在2號站點，約為16.9 cm/s，而夏季近底最大的潮流流速則出現在3號站點，約為16.6 cm/s。

a) 夏季底層

b) 冬季底層

c) 夏季表層

d) 冬季表層圖8 M2分潮表層潮流橢圓的季節變化

M2分潮表層的潮流橢圓季節變化特征表現為：1、5—7號站位表層潮流為旋轉流，2號站位為帶有旋轉性質的往復流，3、4號站位為往復流。各站點表層潮流橢圓旋向與近底潮流橢圓旋向相同，各站冬夏兩季的潮汐橢圓依舊是旋向相反的(圖8)。在潮流流速分布方面，分潮表層潮流隨季節變動不大，冬夏兩季潮流最大流速相近，夏季最大表層流速在2號站點，約為46.5 cm/s；冬季最大表層流速也為2號站點，約在43.5 cm/s。

4.4.2K1分潮

從K1分潮近底層潮汐橢圓對比(圖9)可見，除1號站由夏季往復流變為冬季旋轉流之外，近底層潮流類型幾乎不隨季節發生變化，2—4號站點近底潮流表現為往復流，5號站由于K1分潮強度較弱，較難在圖上表現，其分潮實際上也表現為往復流，而6、7號站點均表現為旋轉流。在潮流橢圓旋向方面，各站夏冬兩季潮流橢圓旋向相反，1—3、5號站的潮流橢圓旋轉方向由夏季的順時針旋轉變為冬季的逆時針旋轉，6、7號站則有夏季的逆時針旋轉變為冬季的順時針旋轉，這與近底層M2分潮的潮流特征相一致，而4號站點的橢圓旋向在夏冬兩季均表現為逆時針旋向。在潮流流速分布方面，1—3號站點近底層冬季與夏季的潮流流速十分接近，4號站點，冬季潮流速度大于夏季，而6—7號站點則表現為夏季大于冬季。珠江河口水域最大近底流速出現在2號站點，冬季約為12.6 cm/s，而夏季約為12.7 cm/s。

從 K1分潮表層的潮流橢圓季節變化(圖9)可以看到，1—5號站的潮流類型不發生變化，其中1、4、5號站表現為旋轉流，2、3號站位為往復流。6、7號站則在夏季表現為往復流，冬季為旋轉流。從橢圓旋向方面看，各站冬夏兩季的潮汐橢圓旋向相反，具體表現為1—5號站夏季為順時針旋轉，冬季為逆時針旋轉；6、7號站則為夏季逆時針旋轉，冬季順時針旋轉。在潮流流速分布方面，1、2、4、6、7號站點表層潮流流速為夏季大于冬季，而3、5號站點表層潮流速則為冬季大于夏季，2號站點有著夏季最大表層潮流流速，約為41.0 cm/s，也有最大的冬季表層潮流流速，約為37.6 cm/s。

a) 夏季底層

b) 冬季底層

d) 冬季表層圖9 K1分潮表層潮流橢圓的季節變化

5 模型應用

隨著海上工程技術的不斷提高，海上橋梁建設和島礁建設等海上工程項目日益增多，然而這些海上工程項目都會對周遭海域的流場產生影響。比如,港珠澳大橋的建設對橋軸線1 km范圍內的潮位，潮流的流態、流速，以及珠江河口的進出潮量都有一定的影響[17-21]。本文在建立河口區數值模型的基礎上進行相關領域的定性數值實驗研究，在模型中添加人工島環境并模擬計算區域的潮流變化。人工島的經緯度為113.65°E，22.25N，面積約為3 km2，人工島的具體位置見圖10。

圖10 模擬人工島位置

為了更加清晰地展現海上工程對于潮流的影響，分別模擬了珠江口洪季和枯季2個時期添加試驗人工島的水動力環境并將其與初始環境相對比。根據珠江口的年徑流量特征，取6—9月份珠江河口洪季，洪季徑流量的平均值為10 600 m3/s，風場數據采用香港機場觀測站6—9月的觀測數據；取12月至次年2月份為枯季，枯季平均徑流量為2 880 m3/s，風場數據采用同一時間段香港機場觀測站的觀測數據。在一個潮周期內，分別對水域水體的表層、中層、底層的潮流場以及大潮期和小潮期的表層潮流場進行分析，并制成平均流速矢量圖，其中，小潮期表層的潮流場與原表層的潮流場變化基本一致，故不作表述。

5.1 表層潮流變化

枯季表層潮流的流速方向為南向，施工后潮流在人工島區域附近發生繞流，人工島東側海域流向發生逆時針偏轉，繞過人工島后發生順時針偏轉；西側發生順時針偏轉，繞過人工島后發生逆時針偏轉，最終流速矢量方向與人工島建設前保持一致(圖11)。取人工島南側A點分析流速，發現其在人工島建設前最大流速為1.26 m/s，在人工島建成后為0.45 m/s，流速減小了64.3%。洪季時表層潮流流速矢量的變化與枯季保持一致。A點的最大流速由建設前的1.30 m/s變為建成后的0.47 m/s，流速減小63.8%。

a) 枯季

b) 洪季注：黑色箭頭表示施工后，淺灰色箭頭表示施工前圖11 施工前后表層潮流平均流速變化

5.2 中層潮流變化

相比于表層較為穩定的潮流場，中層的潮流場發生了劇烈的變化。潮流的流速矢量除了發生繞流現象以外，還在人工島的東南側、人工島的西北側、人工島與西側島嶼間的區域均產生環流現象，使得附近水域的潮流場環境變得更加復雜。例如顯示區域的西南角的流速矢量受到環流影響發生了明顯的逆時針偏轉，而在表層時這部分的流速矢量方向基本不發生變化(圖12a)。環流的特征和潮流場的變化在洪季時表現的更加明顯(圖12b)。此外A點的最大流速也發生了大幅減小的現象，枯季時由建設前的74.82 cm/s變為建設后的30.12 cm/s，流速減小了59.7%；洪季時由建設前的96.52 cm/s變為建設后的40.67 cm/s，流速減小了57.9%。

a) 枯季圖12 施工前后中層潮流平均流速變化

b) 洪季注：黑色箭頭表示施工后，淺灰色箭頭表示施工前續圖12 施工前后中層潮流平均流速變化

5.3 底層潮流變化

底層潮流流速矢量反向與表層發生近180°偏轉，表現為北向。潮流在人工島處發生繞流現象，人工島東側潮流矢量發生順時針偏轉，繞過人工島后發生逆時針偏轉。這與表層的流速矢量變化相似。但西側的潮流流速矢量方向并未發生逆時針偏轉從而繞過人工島，而是產生了順時針方向的偏轉，并在流過人工島區域后又產生順時針的偏轉，這可能是受到了人工島西側產生的環流，以及西側海岸線和島嶼的地形因素的影響(圖13)。此外，A點枯季時的最大流速由建設前的40.20 cm/s減小為13.41 cm/s，減小了66.6%；洪季的最大流速由建設前的41.33 cm/s減小為12.82 cm/s，減小了69.0%。

5.4 大潮期表層潮流變化

大潮期表層潮流場的主要特征為：人工島的西北側會產生新的環流結構，并影響周圍的潮流流動，使得人工島在西側和北側的流速矢量變化異于整個潮流周期內的表層流速變化，這種環流結構及其影響在洪季表現得更為強烈(圖14b)。

a) 枯季

b) 洪季注：黑色箭頭表示施工后，淺灰色箭頭表示施工前圖13 施工前后底層潮流平均流速變化

a) 枯季

b) 洪季注：黑色箭頭表示施工后，淺灰色箭頭表示施工前圖14 施工前后大潮期表層潮流平均流速變化

6 結語

本文采用三維數值模型模擬了珠江河口區域潮流場，通過與觀測數據比對發現模型能夠準確地模擬珠江河口水域的潮流特征。在此基礎上，取特定點進行潮流分析，得到了以下結論。

a) 珠江河口水域的潮型系數約為1.3～1.8，屬于不正規半日潮，M2分潮居于主導地位。珠江河口區域受徑流下泄和狹窄地形的影響，潮流以往復流為主，旋向在近陸架淺海海域主要以旋轉流為主。

b) 河口區潮流場的特征較穩定，季節變化性弱，潮流的性質和類型基本保持一致；僅在潮流橢圓旋向上表現為冬季與夏季相反的特征。

c) 人工島區域附近的流場流速明顯變小，方向發生變化產生繞流現象。人工島與附近島嶼及海岸線之間的相會影響會使區域的流場產生環流現象，進而影響區域的潮流變化。而在寬闊水域區域，潮流方向僅發生繞流現象。