珠江河口磨刀门水道枯季盐水入侵特性分析

方神光

(珠江水利委员会 珠江水利科学研究院, 广东 广州 510611)

珠江河口属典型弱潮型河口, 一定条件下, 容易出现海水与河水之间明显的分层现象[1], 海水在表层淡水下方呈楔状沿河底向上游推进, 形成盐水楔或异重流。如在伶仃洋水域, 来自南海大陆架高盐水体在潮汐作用力下, 主要沿伶仃水道南段和暗士顿水道上溯[2], 与虎门口下泄径流形成混合和分层现象。珠江河口八大口门中, 以磨刀门输水和输沙量最大, 磨刀门水道上游是中山、珠海和澳门的主要水源地, 近年由于气候变化及人类活动等, 河口咸潮入侵频繁, 严重影响沿岸城市的供水安全。

鉴于磨刀门水道的重要性, 有关该水道咸潮入侵的原因和机理研究成果较多, 如刘杰斌等[3]和包芸等[4]分析实测资料显示, 磨刀门水道出现枯季小潮期间咸潮快速上溯, 是由小潮期间磨刀门水道净泄量为负所致; 闻平等[5]认为影响磨刀门水道咸潮入侵的主要因素是径流, 并推荐了最小压咸流量范围以及最佳压咸补淡的时机; 陈荣力等[6]和罗丹[7]分析实测资料, 认为大潮期间加大上游下泄流量是最好的压咸时机; 韩志远等[8]的研究显示, 近些年的河道采砂以及河口围垦导致磨刀门口门“调淡”作用的丧失, 才是咸潮影响加剧的主要原因; 河口咸潮入侵与径流和潮汐作用的此消彼长关系密切, 贾良文等[9]的分析显示, 磨刀门水道枯季潮流为主要动力,河口下层有反向水流, 存在明显的因密度差而形成的密度环流。在磨刀门口, 当上游西江梧州来流量小于2 500 m3/s时, 属于以潮为主河口, 潮流动力作用加强, 咸潮对上游的影响逐步加大。以西蒙斯咸淡水混合判别指标分析, 咸淡水混合状态多为弱混合型及缓混合型, 属高度分层型, 形成明显盐水楔[10]。

可见, 枯季磨刀门水道由于径流减少, 盐水入侵容易形成分层现象, 具有典型的三维水动力特性。鉴于此, 依托磨刀门水道最新地形资料, 本文采用Delft3D软件建立了磨刀门水道三维潮流和盐度数值模型, 并采用2009年12月10~25日的实测资料对模型进行了率定和验证。在此基础上, 对磨刀门水道的潮流和盐度入侵特性进行了分析和探讨。

1 数学模型的建立

1.1 控制方程

控制方程由连续方程、动量方程和盐度方程组成, 在 Delft3D中, 采用曲线贴体σ坐标系, 所建立的数值模型中已考虑盐度正压和斜压对水流运动的影响, 以便更真实地模拟河口径流和潮汐的相互作用。其通用形式如公式(1)。

其中,φ表示流场中不同物理量的通用变量,Rφ(xi)表示源项。uj是贴体坐标系上xj的速度分量;ρ是水体密度, 随温度和盐度变化, 采用国际通用海水密度计算公式(UNESCO), 模型中已考虑由盐度变化引起的水体密度梯度对水流运动的影响;Γφ为扩散项参数。DELFT3D软件中的相关控制方程详细表达形式可参考相关文献[11-12], 并采用交替方法对控制方程进行离散和求解。

1.2 计算区域及网格

图1 磨刀门水道实测站点分布图Fig.1 Sketch of field station locations in Modaomen Waterway

2009年 12月10~25日(农历 10月 24日~11月10日), 珠江水利科学研究院[7]对磨刀门水道竹银-磨刀门口门(大横琴水文站以南约3.9 km)之间近40 km的磨刀门水道中的 8个站点的潮位、流速和盐度进行了逐时测量, 实测点位置如图1所示, 其中1#测点位于磨刀门口门位置, 距大横琴水文站以南约3.8 km,8#站点位于竹银, 4#站点布置在洪湾水道内。数值模拟试验中, 为方便验证和提取边界, 选取的模型计算区域与实测区域一致(图1)。采用正交曲线网格对计算区域进行剖分, 共计 216×27个网格, 垂直方向平均分成10层。

1.3 边界条件及参数选取

1.3.1 边界条件

计算范围内的磨刀门水道边界条件相对较为简单, 以该水道上游竹银位置为上边界, 应用2009年12月10~25日实测流量数据进行控制, 下边界主要是磨刀门河口边界和洪湾水道边界, 采用同一时期实测潮位控制。

1.3.2 初始条件

DELFT3D 软件中, 初始条件有“冷启动”和“热启动”之分, “冷启动”即给定一个零初始场, 依靠模型不断迭代计算以接近实际情况; “热启动”即给定初始水流、水位、盐度等初始场。在潮位和潮流的数值模拟计算中, 边界影响能迅速传递到流场内部,因此模拟的流态能迅速趋向于真实值; 但在盐度场计算中, 由于盐度输移和扩散缓慢以及遵循物质守恒规律, 盐度若给为零初始场, 则要花费相当长的计算周期才能达到真实分布状态。因此, 如果给定的盐度场能尽量接近真实情况, 能提高收敛速度, 节省计算时间。鉴于此, 本文依据2009年枯季磨刀门水道布置的 8个实测点位的盐度资料以及实际地形状况, 采用插值方法给出了计算区域每一层网格节点上的盐度初始值。

1.3.3 其他参数选取

根据模型调试显示, 模型糙率变化范围为0.016~0.025, 从上游至口门依次减少。由于水平方向相对垂直方向大的多, 因此水平方向的紊动黏性系数和盐度扩散系数根据模型调试后的结果, 直接给定为常值即可; 在垂直方向, 则通过水动力k-ε(紊动能-紊动耗散率)紊流模型计算给出。

1.4 数学模型验证

由于选取1#、8#和4#站点分别作为磨刀门口门、磨刀门水道上游以及洪湾水道的边界控制站点, 因此该3个测站不参与验证。图2~图6给出了2#和6#测站的潮位、表底层流速和盐度的验证图, 其他测站由于篇幅所限不一一列出。验证结果显示, 从潮位、流速和盐度各项实测值与数值模拟结果比较来看,数值模拟计算的潮位和流速与实测值吻合较好, 潮位的绝对误差在 0.07 m以内, 其中 2#的误差最小,6#站的误差最大。流速除了个别站位存在较大的误差, 误差均在10%左右。盐度吻合相对较差, 个别实测点底层盐度数值模拟结果与实测误差较大, 尤其是底层靠近上游的验证站位, 这可能是由该测量站位资料存在一定误差引起。但从总体来看, 模拟计算各测点盐度变化总体趋势与实测基本一致, 综合考虑和分析实测过程中各种干扰因素对盐度测量值的影响、局部地形的差异以及数值模拟的误差精度等,本模型基本能反映 2009年 12月 10~25日(农历 10月24日~11月10日)期间的潮位、流速和盐度的变化, 可以应用本模型来分析该水道潮流和盐度变化规律。

图2 潮位验证图Fig.2 Verification of tidal level

图3 2#站流速验证图Fig.3 Verification of velocities in No.2 station

图4 6#站流速验证图Fig.4 Verification of velocities in No.6 station

图5 2#站盐度验证图Fig.5 Verification of salinities in No.2 station

图6 6#站盐度验证图Fig.6 Verification of salinities in No.6 station

2 磨刀门水道潮流特性分析

2.1 流态特性

图7给出了数值模拟的大潮涨落急时刻表层和底层的平面流态(中层及中潮和小潮时刻流态与大潮基本一致, 此处没有给出), 图8和图9给出了实测情况下的表层和底层流速沿河口向上游方向的变化线,横坐标表示以口门 1#测点为起始原点, 向上游方向距离为正。从平面流态来看, 靠近东侧为磨刀门水道主槽, 西侧主要为浅滩分布, 因此东侧流速要明显大于西侧。大潮涨急时刻, 从口外进入磨刀门水道的涨潮流向北, 即西北方向; 洪湾水道水流为东南向,即磨刀门水道部分涨潮流进入到该水道。从表、底分层来看, 表层平均流速略大于底层平均流速。大潮落急时刻, 落潮流速显著大于涨急流速, 磨刀门水道落潮流为东南向; 洪湾水道流向则为西北向, 部分潮流进入磨刀门水道。流速存在分层现象, 由表向底流速依次减小。对中潮和小潮的分析显示, 磨刀门水道表层落潮流速显著大于涨潮流速, 大潮和中潮时, 底层落潮流速略大于涨潮流速, 小潮时, 底层落潮流速略小于涨潮流速; 磨刀门水道表层平均涨潮流速值随潮型增大而增大, 显示涨潮时潮流在磨刀门水道内占优势; 落潮流时, 径流作用凸显, 与落潮流相互作用, 导致河道流速变化较为复杂, 随潮型变化趋势不明显; 各潮型下, 底层涨、落潮流速都随潮型增大而增大。总体来看, 枯季磨刀门水道总体涨、落潮流速都不大, 各潮型下的表层总体涨潮平均流速都在0.5 m/s以内, 总体落潮平均流速在0.8 m/s以内; 底层总体涨落潮平均流速都在0.5 m/s以内。

2.2 表底流态反向现象

枯季小潮期间, 径流量小, 口门潮汐动力弱, 容易出现盐水入侵造成的密度分层现象, 表层淡水向海漂移, 而底层盐水由于补偿流和盐度斜压等作用,向上游运动。图10给出了小潮时刻(12月 11日 01时小潮低潮)磨刀门水道的纵剖面流速等值线, 此时流速处于涨落潮的交替时刻。纵坐标为高程, 坐标原点设为口门位置 1#测点, 横坐标为以口门 1#测点为原点向上游方向的距离。图中流速为正表示流向向海, 负值表示向陆地。可见, 小潮低潮时刻, 底层水流以向陆为主, 表层水流以向海运动为主, 表底层流向相反; 垂向上, 流向变化的位置流速等值线分布较为密集, 流速梯度变化大。实测期间, 该现象在口门 1#测点能明显观察到。大、中、小潮型下, 1#测点表层平均落潮和涨潮历时分别为8.8 h和3.9 h、12.5 h和3.8 h、8.6 h和3.9 h; 底层平均落潮和涨潮历时分别为8 h和16.8 h、7.3 h和17 h、6 h和18.9 h。该测点表层和底层涨落潮时段不完全重合, 说明表层和底层流向存在不一致的时段; 从各潮型来看,表、底层流速相反的时间段一般出现在经过第一个长时间的落潮流之后的第二个较短落潮时间段内。且沿水深方向流向发生变化的位置随潮动力减弱而有向水面移动的趋势。

图7 磨刀门水道大潮期间涨落急流态Fig.7 Flowing at Modaomen waterway during spring

图10 磨刀门水道小潮纵剖面流速图Fig.10 Velocity contour at a longitudinal cross section during neap

3 磨刀门水道盐水入侵特性分析

3.1 盐度平面分布特征

图11给出了数值模式计算得到的大潮涨落急时刻表层和底层的盐度分布, 图12给出了观测站点表层和底层盐度分布。由于磨刀门水道从口门至与洪湾水道汇合口之间河段的主槽靠近东岸, 因此表底层盐度入侵受地形影响明显, 总体呈现: 涨潮时磨刀门水道东侧盐度大于西侧, 落潮时东侧盐度小于西侧的趋势。涨急时, 表层和底层盐度为10的等值线分别到达挂定角和大冲口水闸附近; 落急时, 表层和底层盐度为10的等值线线则都到达灯笼水闸附近。可见, 涨潮时, 盐度为 10的等值线入侵距离超过表层; 落潮时, 两者相差不大; 实测和模拟结果都显示, 磨刀门水道表层盐度显著小于底层盐度,小潮期间该现象最为明显, 且小潮期间的底层盐度要大于大潮和中潮, 显示磨刀门水道枯季小潮期间的盐水入侵更为严重, 并形成明显的盐水分层现象,而表层盐度变化规律仍基本遵循随潮型增大而增大的规律。竹排沙浅滩东侧水道盐度总体小于西侧, 大排沙浅滩东侧石岐水道盐度则明显大于西侧磨刀门水道。从表、底层盐度等值线分布来看, 涨落潮时底层盐度等值线分布梯度显著大于表层, 同一位置的底层盐度明显高于表层。

3.2 盐度垂向分布特征

为探讨大、中、小潮下磨刀门水道盐度在垂向的变化规律, 图13给出沿主槽对应剖面涨落潮下的盐度垂向分布, 纵剖面线位置如图1所示。由于磨刀门水道地形变化对盐水入侵影响明显, 为分析方便,根据水深地形特征, 大致可以将磨刀门水道近口门河段分为4段: 河段1为从口门至与洪湾水道汇合口,河段长度大致为15 km, 从口门向河道上游方向14 km范围内, 河底地形高程从大约–4 m下降至接近–12 m,然后在1 km内迅速上升到约–3 m; 河段2为从洪湾水道汇合口到竹排沙, 该河段长约16 km, 呈现两端高、中间低的形态, 该河段下游端地形高程为–3 m,上游端地形高程接近 0, 中间最低位置地形高程为–12 m左右; 河段3为石岐水道, 该水道位于大排沙的东侧, 河道长度约6 km, 从纵剖面地形上看, 水深较浅, 也呈现两端高、中间低的形式; 河段4为从磨刀门水道进入石岐水道分流口往上游方向河段, 该河段河床地形由–2 m高程在5 km距离内迅速下降到–15 m。根据各潮型下的纵剖面盐度分析各河段盐度特性可见:

图11 磨刀门水道大潮期间盐度分布Fig.11 Salinity contours at Modaomen waterway during spring

1) 大潮涨急时, 外海高盐度(盐度在20以上)的纯海水进入到磨刀门水道河段 1内, 表层和底层入侵距离分别达4 km和10 km左右, 底层高盐水体入侵距离显著大于表层, 高盐海水与冲淡水接触锋面盐度等值线分布密集, 盐度等值线近似垂向分布,与洪湾水道汇合口处盐度约为 12; 河段 2底层盐水入侵距离大于表层, 末端竹排沙位置盐度为 4左右;河段3末端西河水闸位置盐度约为1; 河段4从西河水闸向上游, 盐度都小于1。

2) 大潮落急时, 河段 1表层高盐海水基本退出口门, 为冲淡水占据, 底层由于地形原因, 部分高盐水体聚集在底部无法退出, 形成显著的表底分层现象, 与洪湾水道汇合口处盐度为 16左右, 明显高于涨潮时的盐度; 河段2表、底层盐度向上游入侵的距离反而大于大潮涨急时, 该河段末端竹排沙位置盐度为 6左右; 河段3整体盐度较涨急时有显著升高,末端西河水闸位置盐度约为5; 河段4的整体盐度也呈现显著升高。

图12 磨刀门水道枯季表层和底层涨落潮平均盐度分布Fig.12 Average salinity curves at water surface along Modaomen waterway

图13 不同潮型下的纵剖面盐度等值线Fig.13 Vertical section of salinity contours with different types of tides

3) 中潮涨急时, 河段1底层约7 km范围内由外海高盐度水体占据, 表层则由冲淡水覆盖, 与洪湾水道汇合口位置的最大盐度为5左右; 河段2在灯笼水闸以下游河段的盐度变化范围在 1~5, 底层最远入侵距离达到灯笼水闸, 灯笼水闸往上游河段 2部分、河段3和河段4的水体盐度都小于1。

4) 中潮落急时, 河段 1为冲淡水占据, 呈现显著的表、底分层现象, 底层最高盐度达 17以上, 表层盐度为 4左右, 与洪湾水道汇合口位置的盐度达到6左右; 河段2的盐度较涨急时整体有所升高, 底层 1的等值线向上游最远接近联石湾水闸, 从联石湾水闸往上游河段2的部分、河段3和河段4的水体盐度都小于1。

5) 小潮涨急时, 河段 1底层高盐水体入侵较为显著, 尽管底层含盐度20以上的纯海水仅仅进入到离口门 4 km的范围, 但河段 1的底层基本都被 15以上的高浓度盐水占据, 表层为冲淡水占据, 形成显著分层现象, 与洪湾水道汇合口位置最高盐度为8左右; 河段 2在联石湾水闸下游河道的水体盐度变化范围在 1~10, 较高含盐度出现在洪湾水道汇合口位置底层, 显示是河段 1底层部分较高含盐度的水体侵入所致; 联石湾水闸以上河段的水体盐度都在1以内;

6) 小潮落急时, 河段 1水体含盐度较涨急时有显著下降, 底层最大盐度为 15左右, 表层最低为 3左右, 形成显著表底分层现象, 底层较高浓度的含盐水体由于地形阻碍无法完全退出口门, 盐度等值线呈现由海向陆的倾斜分布, 与洪湾水道交汇处的最大盐度为6左右; 河段2水体含盐度较涨急时有所下降, 联石湾水闸下游河道水体盐度变化范围为1~6, 从联石湾水闸往上游河道含盐度都在1以内。

综上可见, 近河口的磨刀门水道由于其特殊的水深地形和潮流动力特征, 决定了口门盐度入侵的变化特性:

1) 枯季大潮时, 潮动力占绝对主导作用, 外海高盐度的纯海水从口门显著入侵到磨刀门水道, 对磨刀门水道下泄径流形成了显著顶托作用, 致使大潮时的盐度等值线总体呈现垂向分布; 中潮和小潮时, 潮动力减弱, 径流作用凸显, 密度较小的淡水或冲淡水由表层向下游流动, 底层高浓度的盐水由于地形及补偿流动原因向上游入侵, 形成显著的表底分层现象。

2) 从口门至洪湾水道的河段 1, 由于水深地形呈现由海向陆倾斜下降的形态, 导致底层容易为高浓度盐水占据, 且不易完全退出口门; 与洪湾水道的汇合口处地形的抬高对来自外海的底层高盐水体入侵具有显著的阻挡作用。

3) 大潮时的盐水入侵距离较远, 盐度为 1的等值线能达到竹银以上, 中潮和小潮时则基本都在联石湾水闸以下。

4) 大潮和中潮期间, 落潮时的盐水向上游入侵距离反而较涨潮时更远, 根据对流态分析显示, 潮汐动力强时, 涨潮动力强劲且枯季淡水径流动力较弱, 致使上游下泄的淡水径流聚集在磨刀门水道上游河段内无法下泄, 导致上游水道含盐度迅速降低;落潮时, 一方面聚集在水道内的大量冲淡水由表层迅速下泄, 底层高盐水体由于补偿流动向上游入侵;另一方面磨刀门水道水深地形变化迅速, 造成表底层水流的强烈紊动, 进一步加剧了底层高含盐水体往上游方向扩散。小潮期间, 由于整个水道内水流流速很小, 流态平缓, 紊动较弱, 总体仍呈现涨潮时入侵距离大于落潮。因此, 枯季磨刀门水道盐水入侵的主要影响因素为地形和潮动力。

4 结论

为探讨枯季磨刀门水道潮流和盐水入侵特性,本文建立了磨刀门水道三维潮流和盐度数学模型,对2009年枯季近半月的潮流和盐水运动进行了模拟和分析, 结果显示:

1) 枯季由于上游径流量小, 磨刀门水道总体涨、落潮流速都不大, 表层总体涨潮平均流速都在0.5 m/s以内, 总体落潮平均流速在0.8 m/s以内; 底层总体涨落潮平均流速都在0.5 m/s以内;

2) 枯季磨刀门水道表、底层流速反向的时间段一般出现在经过第 1个长时间的落潮流之后的第 2个较短落潮时间段内, 沿水深方向流向发生变化的位置随潮动力减弱而有向水面移动的趋势;

3) 从盐度的平面分布来看, 磨刀门水道近口门河段(与洪湾水道汇合口以下)总体呈现涨潮时, 磨刀门水道东侧含盐度较西侧大, 落潮时东侧含盐度较西侧小的趋势。竹排沙浅滩东侧水道含盐度总体较西侧小, 大排沙浅滩东侧石岐水道含盐度则明显较西侧磨刀门水道大;

4) 潮汐动力较强时(大潮和中潮), 落潮时的盐水向上游的入侵距离反而较涨潮时更远, 潮汐动力弱(小潮)时, 整个水道内水流流速很小, 流态平缓,紊动较弱, 总体仍呈现涨潮时入侵距离大于落潮时的情况。因此枯季磨刀门水道盐水入侵特性取决于水道地形和潮动力。

[1]谭维炎.盐水楔运动规律的研究述评[J].水科学进展,1994, 5(2): 149-159.

[2]徐峰俊, 朱士康, 王华.伶仃洋水动力环境分析及治理策略探讨[J].人民珠江, 2004(1): 11-14, 27.

[3]刘杰斌, 包芸, 黄宇铭.丰、枯水年磨刀门水道盐水上溯运动规律对比[J].力学学报, 2010, 42(6):1098-1103.

[4]包芸, 刘杰斌, 任杰, 等.磨刀门水道盐水强烈上溯规律和动力机制研究[J].中国科学: G辑: 物理学力学天文学, 2009, 39(10): 1527-1534.

[5]闻平, 陈晓宏, 刘斌, 等.磨刀门水道咸潮入侵及其变异分析[J].水文, 2007, 27(3): 65-67.

[6]陈荣力, 刘诚, 高时友.磨刀门水道枯季咸潮上溯规律分析[J].水动力学研究与进展 A辑, 2011, 26(3):312-317.

[7]罗丹.磨刀门咸潮测验的实践[J].水利技术监督, 2011,19(5): 21-23.

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

[9]贾良文, 吴超羽, 任杰, 等.珠江口磨刀门枯季水文特征及河口动力过程[J].水科学进展, 2006, 17(1):82-88.

[10]金腊华, 徐峰俊.河口及近海水质模拟[M].北京: 化学工业出版社, 2007.

[11]陆仁强, 何璐珂.基于Delft3D模型的近海水环境质量数值模拟研究[J].海洋环境科学, 2012, 31(6):877-880.

[12]Lesseer G, Roelvink J A, Van Kester J A T M, et al.Development and validation of a three-dimensional morphological model［J］.Coastal Engineering, 2004,51: 883-915.