磨刀门水道咸潮上溯数值模拟及其分析*

刘祖发，丁 波，关 帅，陈记臣，张泳华，于海霞

(中山大学水资源与环境研究中心∥华南地区水循环与水安全广东省普通高校重点实验室，广东 广州 510275)

磨刀门水道咸潮上溯数值模拟及其分析*

刘祖发，丁 波，关 帅，陈记臣，张泳华，于海霞

(中山大学水资源与环境研究中心∥华南地区水循环与水安全广东省普通高校重点实验室，广东 广州 510275)

根据磨刀门水道2006年11月和2007年1月的实测水位、盐度等资料，利用MIKE3对该地区咸潮上溯进行数值模拟。经率定验证后，水位与盐度模拟结果与实测值契合度较高，两者的拟合效果检测指标均较好。对数值模拟结果进行分析，结果表明：盐度的层化受地形、分流比及潮汐涨落变化等因素的影响，且由分层带来的盐度差造成水体密度差，进而产生垂向环流致使盐度出现异常变化；潮通量与盐通量的变化具有良好的相关性，反映出潮水运动的物质运输载体功能。

数值模拟；咸潮；磨刀门；特征分析

咸潮上溯是河口地区普遍存在的现象之一，尤其在人类活动对自然改造剧烈的区域。近些年来，珠江三角洲地区咸潮上溯呈现频率提高、范围扩大的趋势，其中以磨刀门水道最为严重。磨刀门水道承担着中山、珠海与澳门地区的供水任务，沿岸分布着密集的工农业重镇，经济繁荣且人口密度高，对淡水的需求量很大。连续几年于枯季发生的特大咸潮，给该地区带来巨大的经济损失和社会影响。为应对咸潮带来的影响，国家已经实施珠江压咸补淡应急调水方案，在一定程度上缓解了供水紧张的局面[1-3]。为解决咸潮上溯问题，章文等[4]深入探讨了磨刀门水道盐度变化与潮汐过程的相关性；阮波等[5]精确地模拟了盐水楔的二维流动情况；刘杰斌等[6]研究了磨刀门水道枯季盐水入侵咸界的运动规律。

当海洋大陆架较高盐水团随潮汐涨潮流进入河口，并沿河道向上游推进，盐水扩散、咸淡水混合往内陆上溯，致使内陆河道内水体含盐量超过一定标准即形成咸潮[7]，因而准确把握盐水层化现象及区域通量变化有重要意义。咸潮上溯的研究方法有很多，如实测分析、物理模型[8-9]及经验统计等，目前较广泛应用的为数值模拟[10-13]，如利用FVCOM、POM、EFDC等,但鲜有人采用MIKE3。它是丹麦国家水力研究所开发的水环境管理软件，用于模拟具有自由表面的三维流动系统，包括对流扩散、水质、富营养化和沉积作用等模块。其优点是对三维非恒定流进行模拟的同时，还充分考虑密度变化、水下地形、潮汐变化及气象等条件，能够精确捕捉水流变化的特性。本文利用2006年11月1-15日及2007年1月16-27日天河、灯笼山、联石湾、平岗等站点的实测水位、盐度等资料，对磨刀门水道搭建三维模型进行数值模拟，深入探讨咸潮上溯期间盐度层化影响因素及盐分的输运规律。

1 模型建立

1.1 计算方法

MIKE3的数学基础是雷诺平均化的N-S方程，包括紊流影响以及密度变化，同时包含盐度及温度平衡方程：

(1)

(2)

(3)

(4)

其中，ρ为水的密度；cs为海水中声的传播速度；ui为xi方向的速度分量；Ωij为克氏张量；p为压力；gi为重力矢量；VT为紊动粘性系数；δ为克罗奈克函数；k为紊动动能；T与S分别指温度与盐度；DT与DS分别指相关的温度与盐度扩散系数；t指时间；SS为各自的源汇项。MIKE3水动力学模块采用交替方向隐式迭代法(ADI方法)对质量及动量守恒方程进行积分，且对产生的数学矩阵采用双精度扫描法进行求解[14]。

1.2 模拟范围及网格

本文的建模范围从磨刀门出口伶仃洋20 m等深线处往上游至西海水道天河站，由于区间内小型分汊口的影响可忽略，因此模拟区域为磨刀门水道主干道，包括竹银、平岗、联石湾、灯笼山、挂定角、三灶等测站(见图1)。MIKE3 FM采用非结构化三角形网格，可以在不同地区搭建不同分辨率的网格，通过局部调整的方式协调计算精度与效率，如提高内河道较窄以及河网较为复杂地区的分辨率，或对外海地形较为平坦的地段采用较低分辨率。本次通过试验，可在外海开边界一带采取3 km的分辨率、内河道较宽部分用200～280 m的分辨率，部分细窄河道用70 m甚至40 m分辨率来达成计算精度与计算效率的平衡，共计生成三角网格25 607个，节点数为14 746。

图1 磨刀门水道数值模型计算网格Fig.1 Computational grids of Modaomen waterway numerical model

1.3 边界条件及参数设置

本次模拟下游边界为伶仃洋外海20 m等深线，上游边界为内河道天河站处，垂向采用σ坐标，共有5个计算层。水动力方面，下游边界水位由三灶站实测水位得到，上游边界水位由天河站实测水位给定。盐度场方面，下游边界底层盐度设置为32‰，由于表层海水的盐度受内河道淡水的冲淡影响，因此表层盐度略低于32‰，中间各层在底层与表层间插值输入，上游边界盐度设置为淡水0.001‰；由于资料的有限性，没有考虑风场影响。

1.4 初始条件

初始水位条件对模拟结果的影响有限，可取为模拟范围内各测站实测值的平均。对于盐度初始值，若盐度场模拟达到一定时间，可认为盐度场基本达到稳定，其结果已不受盐度初始条件的影响。在本次模拟中，初始盐度值先设为0，依据区域内某点的盐度过程线来判断盐度场是否趋于稳定，而后在盐度场基本稳定的情况下，取此时的盐度场作为真正的初始盐度条件再次进行计算，得出较为准确的模拟结果。

2 工况选取及模型验证

2.1 工况选取

本次研究共选取2个工况：工况1的时段为2006年11月1-15日，工况2的时段为2007年1月16-27日。本文利用工况1来进行模型参数的率定，然后利用工况2进行模型的验证。

2.2 结果验证

对工况1进行模拟，并适当调整扩散性系数、底部糙率及涡粘性系数等参数，得出与实测值较为相符的模拟结果，其中率定站点为灯笼山、联石湾及平岗站点。利用工况1的参数值，对工况2进行模拟，并对结果进行验证。其中，水位的验证站点选取灯笼山和平岗，盐度验证站点选取灯笼山、平岗和联石湾。由于在实际测量中，盐度测点的位置通常位于水面下约1.5 m处，因此本文采用1.5 m水深所在分层的盐度输出结果进行验证，验证结果如图2和图3所示。

由图2和图3可以看出，各站点的模拟结果与实测值较为吻合，且时间变化基本一致。在此次结果检验中，我们采用均方根误差 (RMSE)及相关性系数来表征拟合程度的优劣，结果如表1所示。水位验证站点的相关性系数分别为0.980 4、0.988 5，两者的RMSE值均极小；盐度验证站点的相关性系数均为0.8左右，考虑到现有的盐度测量会因技术方法问题而存在一定的偏差，盐度模拟结果还是较为理想的。总体来讲，通过检验工况2的模拟结果，文中构建的模型能够较为真实地反映水体的运动特征和盐度变化的基本规律。

图2 磨刀门水道灯笼山站和平岗站水位验证图Fig.2 Water level verification at Denglongshan and Pinggang stations

图3 磨刀门水道灯笼山站、联石湾站和平岗站盐度验证图Fig.3 Salinity verification at Denglongshan, Lianshiwan and Pinggang stations

表1 磨刀门水道数值模拟效果检测值

Table 1 Indicators detection of Modaomen waterway numerical model

模拟变量验证站点RMSE相关性系数r水位/m灯笼山0 11230 9804平岗0 02400 9885盐度/‰灯笼山2 19190 7262平岗0 82280 8386联石湾1 54760 8772

3 磨刀门水道层化特征

层化与混合是河口地区的重要物理过程，控制着水体中物质的垂向分布，也影响着物质输运。磨刀门水道咸淡水的混合，不仅与径流、潮汐因素有关，还受风、波浪、地形等因素的影响，物理机制十分复杂[15-16]，因而呈现出显著地盐度分层现象。根据Prandle[17]的研究，河口的固定点的盐度分层程度可以用下式表示：

(5)

式中，Si为分层系数，δS为表底层盐度差，为垂向平均盐度。

当Si>0.32时，河口是层化状态；当Si<0.15时，河口处于完全混合状态；如果Si介于这两者之间，河口是部分混合状态。在本次研究时段内，根据磨刀门的地形以及外力的作用特征，选取一系列特征点来分析不同因素对磨刀门水道盐度分层状况的影响，各特征点位置如图4所示。

3.1 涨落潮与层化

选取特征点A、B和C，计算其分层系数的逐时变化时间序列，并结合潮位、盐度的变化过程，分析三者的内在联系，图5为特征点C处的潮位、盐度以及分层系数之间的关系图。

图4 磨刀门水道不同断面和点位置图Fig.4 Position of the different sections and points of Modaomen waterway numerical model

如图5所示，磨刀门水道潮汐类型为不规则半日潮，即一天之内潮位有两高两低的不等现象。由于潮汐是盐水入侵的动力源，因而盐度也存在相似的变化，且滞后时间约1～2 h。在日周期内，盐度的最大值出现在高高潮涨憩阶段，最小值由于多种因素影响，出现在高低潮或低高潮的落憩阶段。

分层系数随着潮位的变化出现显著的周期性波动，且其极值出现的时间落后于潮位极值约1～2 h。考虑到潮位与盐度之间也存在类似的关系，我们将盐度和分层系数时间序列进行对比，发现两者的极值出现的时间差基本消失。进一步分析可以发现，无论是大潮还是小潮，其潮水落憩附近皆会出现分层系数的极大值，但高低潮过后出现的分层系数极大值要高于低低潮；在接近大潮涨憩时刻会出现一个较高的极值。对于分层系数的极小值，潮差较小的半日潮一般出现在涨憩时刻附近，对于潮差较大的半日潮，则会在一个周期出现两个极小值，分别处于潮位由低到高和由高到低的变化过程中。磨刀门水道盐水分层的变化，随着潮位变化存在周期性的摆动，不断地在层化与混合之间进行转换，反映了河口区潮汐动力与径流交替混合作用的实质。

图5 潮位、盐度、分层系数关系图(点C)Fig.5 Correlation of tidemark,salinity and stratification(point C)

从空间上进行比较，虽然磨刀门盐度分层随着潮汐的变化都存在周期性的转换，但不同地点的分层与混合程度是不同的，从图6(a)中可以看出：A点的分层系数在低潮时基本都远大于0.32，高潮附近则明显处于混合状态，可认为在整个研究阶段都处于较为明显的层化-混合摆动状态；C点的分层系数虽有微小的波动趋势，但始终小于0.15，可认为在整个模拟期间都处于均匀的混合状态；而B点的变化则介于两者之间，分层系数值的变化范围也很明显，且存在极大值稍大于0.32的时段，而极小值又小于0.15，可认为是在分层与部分混合之间转换，且幅度远远小于A点。从地理位置来看，A点处于磨刀门水道的上游，虽盐度值并不高，但是由于上游处潮汐动力和径流作用此消彼长，交替作用明显，因而存在较为显著的垂向盐度梯度；C点处于咸淡水进行大量交换的口门附近，理论上会经常会出现不同程度的盐水分层现象，但是在我们模拟的大部分时间内，此处潮汐动力的作用十分强烈，而径流作用力相对较弱，整个口门附近皆为较高盐度的咸水，导致其为混合状态良好的咸水；B点处于河道中段，此处河道分叉、浅滩和沙洲较多，极大的影响了潮流的形态，并改变了流速的分布，进而削弱了其分层摆动的幅度。

3.2 河道分汊与层化

选取河道分汊处的特征点D、E，两点所在断面的底程相近，但是D点所在河道比E点宽，因而两河道的分流比不同，造成盐度的分层状态也存在一定的区别，如图6(b)所示。

图6 磨刀门水道分层系数对比Fig.6 Comparisons of stratification indexes at Modaomen waterway

从图6(b)可以看出，在一个潮周期内，涨潮附近两点的分层系数值都小于0.15，可认为都处于良好的混合状态。随着潮周期的发展，两点的盐度分层系数值逐渐上升。其中D点逐渐远大于0.32，因而在整个过程中都处于层化-混合转换状态；而E点的分层系数则相对较小，最大值基本处于0.2附近，趋向于部分混合状态。两点在分层系数上出现非常明显的差异，是由于D点所在河道处流量较大，潮流动力与径流交替作用明显，且影响流速的分层，进而造成盐度分层的差异。

3.3 浅滩与深槽的层化差异

选取分别处于深槽与浅滩位置的特征点F、G，两点的河道高程存在明显的差异，分别约-17、-5 m，导致盐度分层出现明显不同的特征。

如图6(c)所示，和D、E两点的分层差异相似，涨潮附近F、G两点盐度都处于混合状态，但落潮时F点的分层系数远大于G点，存在显著地垂向盐度梯度，而G点仅是部分分层。通过对两点处的流速进行对比，发现深槽的整体流速比浅滩大，尤其是底层流速存在明显的差异。此外浅滩表层落潮流先于深槽出现，底层涨潮流也先于深槽出现。

4 垂向环流

当河道内的水体含有大量盐分将致使其密度发生变化，平均而言比淡水密度大。由于密度较大的盐水由底层入侵河口，可能会引起河流内的垂向环流，即水流表层与底层之间的流速方向相反，表层为下泄流，底层为上溯流。

根据分析结果，在磨刀门水道上选取特征点H。盐度变化受流速直接影响，当流速为正时盐度应上升；流速为负时盐度应下降，两者存在着明显的相关关系。但在对H点处盐度与水位、流速的分析中发现有部分时段例外，如工况2中1月18日13：00-15：00、22日17：00-19：00等，这些时段都处于低低潮—高高潮、高低潮—低高潮转变的时刻，即为盐度迅速上升的时刻。共有特点为速度上显示为落潮，而盐度不但没有下降，反而有所上升。这种异于平常的表现，说明可能存在有更加内在缘由。因此，我们以较为显著的1月18日13：00-15：00时间段分析，其中流向取与河道横断面的夹角，逆时针为负，顺时针为正，如图7所示。

从图7可看出，在13：00-15：00这段时间内，断面流速为负指向下游，理论上表层盐度应该的不断下降，但事实上却处于上升状态。对比表层与底层流场图(见图8)，在这段时间内出现了显著的表底层流向相反的现象(见图9)。考虑到在河口地区转潮因素也会引起表底层流态的差异，但通常只维持20～30 min，且较为紊乱。此外，分析H点处分层状况，在此时间段内其分层系数逐渐达到极大值，由此可判定正是垂向环流的出现才造成该时段内在落潮的情况下盐度反而上升的异常现象。垂向环流对盐度的变化趋势与积累有明显的影响。对照H的水位变化可发现，1月18日13：00-15：00正好处于低高潮至低低潮的转变阶段，长时间的垂向环流加强了表底层盐分的流动，使得表层盐度的降低趋势得以削弱，从而有助于盐度的累计并提高盐度峰值，增强了盐水上溯的幅度。

图7 点H表层盐度异常现象Fig.7 Abnormal changes of salinity at H point in Modaomen waterway

图8 点H处表底层垂向环流流场分布Fig.8 Flow field of surface and bottom layers at Hpoint in Modaomen waterway

图9 点H处表、底层流速方向Fig.9 Flow direction of surface and bottom layers at H point point in Modaomen waterway

5 潮通量与盐通量

准确描述与预测盐分在河口中的输移过程及其通量是解决咸潮上溯问题的基础，对于“通量”的概念，沈焕庭[18]将其定义为：在一定时间内通过某一面积的某种物质的质量或者体积。参考程香菊等[19]的研究方法，若河口区域输运的物质在单位体积水体里的含量为C，则该物质输运通量为：

(6)

式中，x、y分别为断面的切向和法向坐标，x1、x2分别为断面的起点和终点；t为时间；v为断面流速；T为潮周期，此次计算中近似取25 h；D为水深；表示分层数，在文中取为5。上式用于计算潮通量时C=1；计算盐通量时C=ρS，其中ρ为海水密度。

研究工况共包括11个潮周期，其中第4个潮周期(农历腊月初二)为大潮，第10个周期(农历腊月初九)为小潮。在磨刀门水道上选取断面1、断面2和断面3，分析各周期内潮通量与盐通量的变化特征，如图10-11所示。

图10 潮通量与盐通量关系图(断面2)Fig.10 Correlation of tidal flux and salt flux (section 2)

潮水运动是河口物质运输的主要载体，因而潮通量与盐通量具有良好的相关性，如图10所示，经过统计相关分析，断面2两者的相关系数为0.85，故其变化具有相似的规律两者都在小潮过后表现为向上游流动，大潮过后变为向下游运动[20]，反映出咸潮上溯的本质为盐分随潮水的迁移。

图11 磨刀门水道特征断面潮通量和盐通量对比图Fig.11 Tidal flux and salt flux comparison of Modaomen waterway sections

从图11(a)可看出，三个断面潮通量变化趋势基本一致。但断面2的潮通量绝对值略高于断面3，原因可能是由于断面B下游存在洪湾水道、鹤洲水道等，对水流有分流作用。此外，在研究的11个周期内，口门内下泄的水量明显远高于上涨的水量，这与实际情况是相符的。

从半月周期角度分析盐通量的变化(见图11(b))，第1、2周处于小潮过后，盐通量迅速上涨且皆为正值；在接近大潮的第3个周期时盐通量变为负值，之后上下波动但始终保持在0刻度线之下，直至第10个周期小潮再次出现，才逐渐越过0刻度线并保持上升。整体来讲在小潮之后盐通量值大于零，表现为上溯量；而在大潮过后小于零，盐水不断下落。从空间角度分析，选取的特征断面中，断面3的盐通量绝对值始终最大，断面2次之，而断面1的始终为最小。参照3个断面所处的地理位置，说明上游的盐通量明显小于下游。故在通过断面的流量相差不大的情况下，下游的水体中含盐量较高，因而盐通量也相对较大。

6 结 论

1)MIKE3适用于磨刀门水道咸潮上溯的三维数值模拟。利用均方根误差、相关系数等进行拟合检验，在水动力方面，水位模拟值与实测值相关系数均在0.98以上；在盐度拟合方面，模拟结果与实测值的契合度很高，真实地反映了咸潮上溯过程流场的性质与盐度的变化规律；

2)磨刀门水道盐度的层化情况受诸多因素影响。固定点的分层系数随潮位的变化出现周期性的层化-混合交替转换，落潮时的分层现象比涨潮时更加明显，且潮汐动力和径流作用交替作用明显的河段断面分层更显著；在地形复杂的地带，水流运动受地形急剧变化的物理因素干扰造成盐度运动复杂易趋向于混合；盐度分层还受河道分汊处的分流比、河道底岸高程等因素影响；

3)层化现象带来的盐度差造成水体密度差，进而在特殊情况下产生的垂向环流致使表层盐度在落潮时出现异常变化，且对盐度的累积上涨起到了不可忽视的作用；

4)潮通量与盐通量的变化具有良好的相关性，反映出潮水运动的物质运输载体功能。

[1] 刘斌，刘丽诗，吴炜，等．基于取淡与流量控制的压咸调度方案研究[J]．水文，2013(4)：84-86;74．

[2] 闻平，戚志明，刘斌．西北江三角洲压咸流量初步研究[J]．水文，2009(9)：74-75;52．

[3] 尹小玲，张红武，方红卫．枯季磨刀门水道咸潮活动与压咸控制分析[J]．水动力学研究与进展A辑，2009(5)：554-559．

[4] 章文，刘丙军，陈晓宏，等．珠江口磨刀门水道盐度变化与潮汐过程的相关性分析[J]．中山大学学报：自然科学版，2013，52(6)：11-16．

[5] 阮波，包芸．磨刀门水道垂向盐水楔及其流动的二维精细模拟[J]．计算机辅助工程，2014，23(4)：35-39．

[6] 刘杰斌，包芸．磨刀门水道枯季盐水入侵咸界运动规律研究[J]．中山大学学报：自然科学版，2008，47(增刊2)：122-125．

[7] 包芸，刘杰斌，任杰，等．磨刀门水道盐水强烈上溯规律和动力机制研究[J]．中国科学，2009(9)：1527-1534．

[8] 方立刚，李宏丽，陈水森．基于反射光谱的珠江口咸潮遥感方法[J]．水科学进展，2012，23(3)：403-408．

[9] 卢陈，袁丽蓉，高时友，等．潮汐强度与咸潮上溯距离试验[J]．水科学进展，2013，24(2)：251-257．

[10] PRITCHARD D W．Salinity distribution and circulation in the Chesapeake Bay estuarine system[J]．Journal of Marine Research，1952，11：106-123．

[11] BOWEN M M．Salt transport and the time dependent salt balance of a partially stratified estuary[J]．Journal of Geophysical Research，2003，108：33-42．

[12] 陈水森，方立刚，李宏丽，等．珠江口咸潮入侵分析与经验模型—以磨刀门水道为例[J]．水科学进展，2007，18(5)：751-755．

[13] PRANDLE D．Saline intrusion in partially mixed estuaries [J]．Estuarine Coastal and Shelf Science，2004，59(3)：385-397．

[14] 马腾，刘文洪，宋策，等．基于MIKE3的水库水温结构模拟研究[J]．电网与清洁能源，2009，2(2)：68-71．

[15] 韩志远，田向平，刘峰．珠江磨刀门水道咸潮上溯加剧的原因[J]．海洋学研究，2010，28(2)：52-59．

[16] 诸裕良，闫晓璐，林晓瑜．珠江口盐水入侵预测模式研究[J]．水利学报，2013(9)：1009-1014．

[17] PRANDLE D．On salinity regimes and the vertical structure of residual flows in narrow tidal estuaries[J]．Estuarine Coastal and Shelf Science，1985，20：615-635．

[18] 沈焕庭．长江河口物质通量[M]．北京：海洋出版社，2001．

[19] 程香菊，詹威，郭振仁，等．珠江西四口门盐水入侵数值模拟及分析[J]．水利学报，2012，43(5)：554-563．

[20] WANG B，ZHU J R，WU H，et al．Dynamics of saltwater intrusion in the modaomen waterway of the Pearl River estuary[J]．Science China(Earth Sciences)，2012，55(11)：1901-1918．

Numerical simulation and analysis of saltwater intrusion to Modaomen waterway

LIUZufa，DINGBo，GUANShuai，CHENJichen，ZHANGYonghua，YUHaixia

(Center for Water Resources and Environment∥Key Laboratory of Water Cycle and Water Security in Southern China of Guangdong Higher Education Institutes，Guangzhou 510275, China)

Saline water intrusions at Modaomen waterway were numerically simulated by using MIKE3 based on measured water levels and salinities in November 2006 and January 2007. Calibrated and validated model showed good fitness for both water level and salinity in case of prediction and observation with good detection indexes. The results indicated that salinity stratification, influenced by topography, diversion ratio and other factors such as tidal fluctuation change, resulted in water density difference and then produced extraordinary change in salinity vertical circulation. The change of salt flux is relative to the tidal flux, suggesting that tide is the carrier of brine.

MIKE3；saline water intrusion；Modaomen waterway；characteristic analysis

10.13471/j.cnki.acta.snus.2016.06.001

2015-12-11基金项目：国家自然科学基金资助项目(41301627)

刘祖发(1961年生)，男；研究方向：水文与水资源研究；E-mail:eeslzf@mail.sysu.edu.cn

P731.23

A

0529-6579(2016)06-0001-09