王光宇，宏波

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

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

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

1 模型介紹

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

動量方程：

?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)

連續(xù)方程：

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

(4)

(5)

狀態(tài)方程：

ρ=ρ(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——垂向紊動黏滯系數(shù)；Kv——垂向紊動擴散系數(shù)；f——科里奧利系數(shù)；ρ——混合密度；ρ0——參考密度；H——總水深；h——未擾動的z坐標原點以下的水深；p——壓力；Sa——鹽度；T——溫度；ξ——自由地勢能；Qu、Qv——動量在x和y方向的源匯項；QT——溫度的源匯項。

2 模型建立

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

圖1 模型計算區(qū)域及網(wǎng)格劃分

3 模型驗證

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

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

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

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

a) 燈籠山

b) 黃埔

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

d) 內(nèi)伶仃

e) 桂山島續(xù)圖3 燈籠山、黃埔、九州港、內(nèi)伶仃、桂山島站水位驗證結(jié)果

a) 三灶

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

a) 水深為4 m的流速

b) 水深為5 m的流速

圖5MO1站洪季流速矢量驗證結(jié)果

a) 水深為4 m的流速

b) 水深為5 m的流速

圖6MO1站枯季流速矢量驗證結(jié)果

4 潮流分析

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

4.1 潮流性質(zhì)

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

表1 各站點表層與底層潮型數(shù)

4.2 潮流類型

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

表2 各站點M2和K1分潮表層與底層潮型數(shù)

4.3 潮流流速的季節(jié)變化

分別統(tǒng)計珠江河口水域底層和表層冬(1—2月)夏(7—8月)兩季的平均潮流流速。從底層的潮流平均流速季節(jié)變化(圖7a)中可以看出，除5號站點外，模型區(qū)域冬夏兩季潮流的平均流速的大小基本保持一致，不發(fā)生明顯變化，流場性質(zhì)較為穩(wěn)定，潮流方向均以北向為主。冬夏兩季最大的平均流速均發(fā)生在2號站點。對比與近底層的潮流流速，模型區(qū)域表層的流速存在一定的季節(jié)變化特性，這與徑流和季風的季節(jié)變化有關(guān)。從夏季到冬季，1、6號站流速明顯變小，3、4、5號站的流速變大，2、7號站點流速大小較為穩(wěn)定。模型區(qū)域的流向發(fā)生明顯的西向偏轉(zhuǎn)，由夏季的偏東方向變?yōu)槎镜钠鞣较颉Ｆ渲?、7號站點的流速方向發(fā)生180°偏轉(zhuǎn)，由東南向轉(zhuǎn)為西北向(圖7b)。相較于穩(wěn)定的底層潮流流速，表層潮流流速的季節(jié)性變化應當與珠江口徑流下泄、季風風向以及外海潮波的季節(jié)性變動有關(guān)。

a) 底層

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

4.4 潮流橢圓的季節(jié)變化

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

4.4.1M2分潮

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

a) 夏季底層

b) 冬季底層

c) 夏季表層

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

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

4.4.2K1分潮

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

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

a) 夏季底層

b) 冬季底層

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

5 模型應用

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

圖10 模擬人工島位置

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

5.1 表層潮流變化

枯季表層潮流的流速方向為南向，施工后潮流在人工島區(qū)域附近發(fā)生繞流，人工島東側(cè)海域流向發(fā)生逆時針偏轉(zhuǎn)，繞過人工島后發(fā)生順時針偏轉(zhuǎn)；西側(cè)發(fā)生順時針偏轉(zhuǎn)，繞過人工島后發(fā)生逆時針偏轉(zhuǎn)，最終流速矢量方向與人工島建設(shè)前保持一致(圖11)。取人工島南側(cè)A點分析流速，發(fā)現(xiàn)其在人工島建設(shè)前最大流速為1.26 m/s，在人工島建成后為0.45 m/s，流速減小了64.3%。洪季時表層潮流流速矢量的變化與枯季保持一致。A點的最大流速由建設(shè)前的1.30 m/s變?yōu)榻ǔ珊蟮?.47 m/s，流速減小63.8%。

a) 枯季

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

5.2 中層潮流變化

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

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

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

5.3 底層潮流變化

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

5.4 大潮期表層潮流變化

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

a) 枯季

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

a) 枯季

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

6 結(jié)語

本文采用三維數(shù)值模型模擬了珠江河口區(qū)域潮流場，通過與觀測數(shù)據(jù)比對發(fā)現(xiàn)模型能夠準確地模擬珠江河口水域的潮流特征。在此基礎(chǔ)上，取特定點進行潮流分析，得到了以下結(jié)論。

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

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

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