珠江河口潮流的数值计算研究
2019-10-17王光宇宏波
王光宇，宏波
(华南理工大学土木与交通学院，广东广州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.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) 人工岛区域附近的流场流速明显变小，方向发生变化产生绕流现象。人工岛与附近岛屿及海岸线之间的相会影响会使区域的流场产生环流现象，进而影响区域的潮流变化。而在宽阔水域区域，潮流方向仅发生绕流现象。
