磨刀门盐水入侵特性及其动力机制研究

吕紫君，冯佳佳，孔 俊，罗照阳，潘明婕

(1.河海大学 江苏省海岸海洋资源开发与环境安全重点试验室，南京 210098;2.南京市水利规划设计院股份有限公司，南京 210006)

盐水入侵是沿海河口地区特有的一种水文现象，多发于枯水期，是河口研究中关键问题之一。国外关于盐水入侵的研究从20世纪50年代开始，早期研究工作主要集中在河口环流、盐水入侵长度、盐淡水混合类型、最大浑浊带等[1-3]。河口区的盐淡水混合与输运是一个极其复杂的过程，盐水入侵受向陆和向海两种盐度输送驱动力共同控制；Hansen[2]较早开始采用机制分解方法研究河口盐度输运，提出以河口流速及盐度的垂、横向变化来描述河口的盐度通量；Lerczak[4]根据实测资料研究了哈德逊河口的盐度输运，将纵向盐度输运分成平流输运、稳定剪切输运及潮汐震动输运。国内关于盐水入侵的研究从20世纪80年代开始，主要集中在长江口[5]和珠江口[6]。

珠江口磨刀门水道沿途分布着众多取水口，是珠海、中山、澳门等城市的主要饮用水水源地。21世纪以来，受航道疏浚、口门挖沙以及上游径流量变化的影响，磨刀门盐水入侵日趋严重，众多学者针对这一问题开展了大量研究。贾良文等[7]对磨刀门水流及盐度的变化过程进行分析，得出磨刀门枯季盐淡水类型属于缓混合型；韩志远等[8]基于不同年份的潮位和地形资料，分析了磨刀门盐水入侵加剧的原因；卢陈等[9]采用物理模型实验方法对不同潮差驱动下盐水入侵距离进行研究，结果表明存在潮差临界值使得盐水入侵距离最短。由于实地测量的时空覆盖有限，而物理试验的研究成本高，所以数值模拟已广泛应用于河口盐水入侵的研究，Gong和Shen[10]采用ELCIRC模型对磨刀门盐水入侵与径流和潮差的响应法则进行研究，表明径流量的大小是影响盐水入侵的关键因素；Wang等[11]、陈文龙等[12]利用FVCOM模型对磨刀门水道的盐水入侵规律进行了分析，发现盐水入侵最远出现在小潮之后的中潮。

上述成果多为对磨刀门盐水入侵规律与影响因素的研究，相应的物质输运通量及其分解的研究偏少，盐水入侵半月周期变化的动力学机制尚不明确，同时磨刀门盐度数值模拟的精度及效率有待进一步提高。鉴于此，本文基于一种最新发展的河口海洋模型SCHISM[13]，构建磨刀门水道三维水流盐度数值模型，利用该模型对磨刀门盐水入侵特性、纵向盐度输运进行分析，探究磨刀门盐水入侵的动力机制。

1 数值模型建立

1.1 数值模型

研究中采用的模型为SCHISM[13]，该模型是基于SELFE[14]开发出的一个跨尺度湖泊-河流-河口-海洋数值模型，计算区域水平方向上采用三角形和四边形混合网格，可以灵活地拟合复杂的地形和岸线；垂向上采用传统的混合SZ坐标或者新增的LSC2坐标[15]，混合SZ坐标不仅能适应复杂不规则地形，而且能有效改善模型垂向的空间分辨率，克服河口区因地形变化大而引起的模拟水平梯度失真的难题；而LSC2坐标在保证深海区垂向网格精度的同时，避免浅水区垂向层数过多。模型采用高阶ELM方法处理动量方程中的对流项，采用迎风格式或隐式二阶TVD格式处理物质输运方程中的对流项，既可提高计算的时间步长，又可保证模型计算稳定性和效率。采用GLS模型求解掺混系数，并提供了与通用的GOTM模块的接口。关于SCHISM模型的其他详细介绍，可见参考文献[13]，模型在水动力模型基础上耦合了波浪、生态、水质、风暴潮等模块，已经用于全球的水动力相关问题的研究与预报，如北海和波罗海[16]等。

1.2 计算区域与网格

磨刀门三维数值模型共有2个开边界，上游开边界位于百顷头以上的外海大桥附近，给流量过程，下游至磨刀门口外-15 m等深线附近，给水位过程；边界流量、水位及盐度均从已建的珠江口二维模型中提取，稳定的初始盐度场通过循环运行模型得出。模型水平方向采用非结构三角形网格，网格共23 283个，节点共13 168个；模型垂向采用混合SZ坐标，平均海平面45 m以下采用Z坐标，分1层；45 m以上采用S坐标，分11层，每层厚度随地形变化。网格分辨率的大小直接影响着模型的计算效率和精度，故分辨率分区域给定，口门区域分辨率取为1 500 m，而为确保盐度模拟精度，河道内分辨率取在40 m以下。

1.3 模型计算参数

模型计算时间为2009年1月2日～1月31日，共30 d，时间步长为150 s；模型全域给恒定曼宁系数，大小取0.015；模型为斜压模型，采用冷启动，不考虑底床变形和风应力，采用GLS湍流混合模型求解掺混系数，采用隐式TVD2格式处理输运方程中的对流项，采用高阶ELM方法处理动量方程中的对流项，干湿网格判断最小水深均为0.1 m。

图1 磨刀门水深地形及测点位置图Fig.1 Topography of the Modaomen estuary and the locations of measurement stations

2 模型验证

采用2009年1月的实测资料对建立的数值模型进行验证，验证点位置分布见图1，其中三灶站、灯笼山站、竹银站为水位测站，M1 ～ M8点为流速和盐度测点。为定量评估验证精度，引入均方根误差(RMSE)，其表达方式如下

(1)

式中：Xmod为模型计算值；Xobs为实际测量值；N为观测次数。鉴于篇幅所限，本文只列出了部分验证成果，如图2～图4所示，圆点为实测值，实线为模拟值。验证表明，计算值与实测值吻合较好，误差均在允许范围之内。总体上，模型能较好地反演磨刀门水流及盐度的变化过程，因此建立的三维水流盐度数值模型可用于盐水入侵机制的研究。

图3 M2、M6测点流速验证Fig.3 Calibration of flow velocity at M2 and M6 station

图4 M2、M6测点盐度验证Fig.4 Calibration of salinity at M2 and M6 stations

3 结果分析

3.1 磨刀门盐水入侵特征

图5 M6点表、底层盐度时间过程线Fig.5 Temporal variation of surface and bottom salinity at M6 station

珠江河口的潮汐类型为不规则半日潮，一天中出现2次高潮和2次低潮，潮差和潮历时明显不等。图5给出了M6点表、底层的盐度时间变化过程，为说明盐度随大小潮潮型的变化，同时给出了三灶站的潮位过程。由图可知，M6点的底层盐度变化在0～9 ppt之间，表层盐度变化在0～6 ppt之间；且盐水入侵过程有着明显的半月周期变化规律，表底层盐度在小潮后期骤然升高，并在中潮时达到峰值，之后随着潮差增大，盐度慢慢变小。

图6和图7分别给出了小潮、大潮期间磨刀门水道潮平均的表、底层的流场和盐度场分布情况，小潮期间(图6)，表层余流在整个磨刀门水道包括洪湾水道及鹤州水道都是向海运动，磨刀门水道流速有明显的横向差异，东侧流速明显较西侧大，这是因为东侧为深槽，是重要潮流通道，流速相对较大，而西侧为浅滩，底摩阻相对较大，削弱了潮流动力，流速相对较小；底层余流在磨刀门水道主要向陆运动，在西岸浅滩区域则向海运动，在洪湾水道上段向海运动，下段则向陆运动；盐水沿磨刀门水道东侧深槽入侵，底部10 ppt盐水可入侵到挂定角以上4 km附近，表层10 ppt等盐度线则处于口门大横琴附近，表底层盐度差异大，盐水分层明显。大潮期间(图7)，表层余流在磨刀门口外海区、口内河道和洪湾水道仍是向海运动，底层余流在磨刀门水道深槽区向陆运动，在浅滩则向海运动；与小潮相比，磨刀门水道内盐水入侵减弱，底部10 ppt盐水回落到挂定角以上1.5 km处，盐水分层较弱。

从潮周期平均的流场和盐度场分布可以看出，磨刀门水道东侧深槽是高盐水入侵磨刀门的重要通道，小潮期间磨刀门水道拦门沙至挂定角段底部积聚较高的盐水，随着潮差增大，积聚的盐水被输送至上游。

图6 小潮期间表层(左)及底层(右)流场、盐度分布Fig.6Surface(left)andbottom(right)meancurrentandtidally-averagedsalinityatneaptide图7 大潮期间表层(左)及底层(右)流场、盐度分布Fig.7Surface(left)andbottom(right)meancurrentandtidally-averagedsalinityatspringtide

3.2 磨刀门盐水入侵动力机制分析

3.2.1 盐淡水混合特征

为了进一步分析磨刀门盐水入侵的动力机制，利用数值模型计算结果，绘制了涨急、落急时刻A-A纵断面的瞬时盐度场，纵断面位置如图1所示。选取分层系数(n)来描述不同时刻盐度的分层特征，n为测点表、底层的盐度差与测点垂线平均盐度的比值，定义为

(2)

图8-a为小潮期间涨急、落急时刻磨刀门纵断面盐度分布图，图上横轴为与外海A点的距离，流速箭头代表沿河道的流速(下同)。小潮期间，三灶站的涨潮潮差为0.93 m，落潮潮差为0.68 m，潮汐动力较小，纵断面的盐度等值线走势比较平缓，沿程盐度分层明显。涨急时刻，口外底部水体盐度变大，最高盐度可达30 ppt，而表层盐度值较低，随着涨潮动力增强，口外高盐水团沿河道底部向上游推进，表层水体受其挤压作用，形成明显的高、低盐水分界面；此时，n2= 0.77，n4= 1.65，n6= 0.01，磨刀门水道的盐水几乎处于高度分层状态，且越往上游，分层越明显，直到上游灯笼山附近，盐水浓度很低，盐度均匀分布，分层现象消失。落急时刻，中层及表层盐水随落潮流退至外海，但口门至洪湾水道段底部仍积蓄着高盐水；此时，n2= 0.47，n4= 0.91，n6= 1.28，水道下游的盐水处于缓混合状态，而洪湾水道至灯笼山的盐水仍处于高度分层状态。

中潮期间(图8-b)，三灶站的涨潮潮差为1.01 m，落潮潮差为1.73 m，与小潮期间相比，纵断面的盐度等值线趋势有所倾斜，盐水分层减弱。涨急时刻，口外盐水楔进入河道，与表层低盐水混合后整体向上游入侵，水道内盐水浓度明显变高，盐水入侵距离增加；此时，n2= 0.51，n4= 0.98，n6= 1.46，在盐水楔锋面所在范围内，盐水处于高度分层状态。落急时刻，口外高盐水稍有退却，盐水入侵距离达到峰值；此时，n2= 0.34，n4= 0.55，n6= 0.35，盐淡水混合充分，河道内的盐水基本处于缓混合状态。

大潮期间(图8-c)，三灶站的涨潮潮差为1.31 m，落潮潮差为2.01 m。涨急时刻，口外高盐水团进入河道，相比中潮涨急时刻，磨刀门水道内盐水浓度变低，盐水入侵距离减小；此时，n2= 0.62，n4= 1.24，n6= 0.70，磨刀门水道的盐水分层加强。落急时刻，口外高盐水随落潮流退出河道，此时，n2= 0.41，n4= 0.89，n6= 0.67，盐水基本处于缓混合状态，大潮潮汐动力强，一天中盐水浓度变化剧烈，盐水“大进大出”。

3.2.2 盐通量分解

本文采用Lerczak等[4]提出的计算盐度通量的方法，将河口的盐度输运分为平流输运(advection transport)、稳定剪切输运(steady shear transport)、潮汐震荡输运(tidal oscillatory transport)三部分。断面总盐通量可以表示如下

FS==<(u0+uE+uT)(S0+SE+ST)dA>

(3)

式中：u为轴向流速，m/s；s为盐度，ppt；A为横断面面积，m2；下标0为潮平均的余流；下标E为剪切流；下标T为潮流波动。MacCready[17]认为公式中的交叉项可以近似为零，得出盐度通量的分解表达式如下

FS=<(u0s0uESE+uTST)dA>=Qfs0+FE+FT

(4)

式中：FS为总盐通量；Qfs0为平流输运项；FE为稳定剪切项；FT为潮汐震荡项。平流输运项是由上游流量以及外海stokes漂流引起的，主要是将上游盐水输运至外海；稳定剪切项主要由垂向剪切流引起，而潮汐震荡项是由流速与盐度的潮内变化相互作用产生，它们则主要是将外海的高盐水输运至河口上游。

以口门处的M2点为起始点，向上游河道每隔5 km取一个横断面，共选取6个横断面，计算并分析各断面的盐通量在一个潮周期过程中的变化情况。图9给出了断面2及断面5的盐通量分量Qfs0、FE和FT以及总盐通量FS在一个完整潮周期过程中的变化情况，图中正值表示向海输运；负值表示向陆输运，断面具体位置如图8所示。

图9-c为断面2的盐通量分量，在断面2处，Qfs0始终为正值，即平流输送作用向海输送盐水，且在小潮期间平流输送作用弱，在大潮期间平流输送作用强；FE始终为负值，即稳定剪切作用向陆输送盐水，且在小潮期间稳定剪切作用强，在大潮期间稳定剪切作用明显减弱；潮汐震荡作用FT在断面2处向陆输送盐水，且在大潮期间最强，在小潮期间最弱；总盐通量FS在小潮期间通过断面2向河道内输送盐水，而在大潮期间向外海输送盐水。

图9-d为断面5的盐通量分量，在断面5处，稳定剪切输送、平流输送和潮汐震荡的规律与断面2基本相同，但输送强度变小；FS在后半月周期的小潮期间向陆输送盐水，随后在中潮转大潮、大潮期间向海输送盐水，而在前半月周期(1月6日～1月17日)，FS基本为0，表明在Qfs0、FE和FT三者作用下，此期间磨刀门水道内的日潮平均总盐量处于平衡状态。

9-a 三灶站水位 9-b 盐水入侵长度

9-c 断面2的盐通量分量 9-d 断面5的盐通量分量图9 盐通量分量过程图Fig.9 Decomposition of salt transport flux

4 结论

基于SCHISM模型，对磨刀门水道盐水入侵特性、纵向盐度输运进行了分析，得出以下结论：

(1)枯季磨刀门盐水入侵一般特性为：小潮期间，磨刀门水道盐水分层明显，拦门沙至挂定角段底部积聚高盐水，盐水入侵距离较小；中潮期间，盐水混合充分，小潮时积聚的底部高盐水逐渐掺混至表层，并被输送至上游，盐水入侵距离最远；大潮期间，一天中盐水浓度变化剧烈，盐水“大进大出”，盐水入侵距离比中潮期间有所减小，但大于小潮期间的入侵距离。

(2)磨刀门日平均盐度受平流输送和稳定剪切输送共同作用，盐量输运在大部分时间内处于不平衡状态。小潮期间，总盐通量向陆，稳定剪切作用强，平流输运作用弱，盐量在河道里累积；中潮期间，总盐通量由向陆转为向海，河道内的日平均盐度最大，盐水入侵距离最远；大潮期间，总盐度通量向海，驱动盐水向陆输送的稳定剪切和潮汐震荡作用小于驱动盐水向海输送的平流输送作用，盐量被冲出河口。

参考文献：

[1]PRITCHARD D W. Salinity distribution and circulation in the Chesapeake Bay estuarine system [J]. Journal of Marine Research, 1952, 11(2): 106-123.

[2]HANSEN D V, RATTRAY M. Gravitational circulation in straits and estuaries [J].J.mar.res,1965, 11(3): 319-326.

[3]SIMMONS H B, BROWN F R. Salinity effects on estuarine hydraulics and sedimentation [C]//Tanaka H. Process of 13th IAHR (3).Kyoto, Japan: International Association for HydroEnvironment Engineering and Research, 1969:192-216.

[4]LERCZAK J A, GEYER W R, CHANT R J. Mechanisms Driving the Time-Dependent Salt Flux in a Partially Stratified Estuary [J]. Journal of Physical Oceanography, 2006, 36(12): 2 296-2 311.

[5]CHEN W, CHEN K, KUANG C, et al. Influence of sea level rise on saline water intrusion in the Yangtze River Estuary, China [J]. Applied Ocean Research, 2016, 54:12-25.

[6]欧素英. 珠江三角洲咸潮活动的空间差异性分析 [J]. 地理科学, 2009, 29(1): 89-92.

OU S Y. Spatial difference about activity of saline water intrusion in the Pearl River Delta [J]. Scientia Geographica Sinica, 2009,29(1): 89-92.

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

JIA L W, WU C Y, REN J, et al. Hydrologic characteristics and estuarine dynamic process during the dry season in Modaomen Estuary of the Pearl River [J]. Advances in water science, 2006, 17(1): 82-88.

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

HAN Z Y, TIAN X P, LIU F. Study on the causes of intensified saline water intrusion into Modaomen estuary of Pearl River in recent years [J]. Journal of Marine Science, 2010, 28(2): 52-59.

[9]GONG W, SHEN J. The response of salt intrusion to changes in river discharge and tidal mixing during the dry season in the Modaomen Estuary, China [J]. Continental Shelf Research, 2011, 31(7): 769-788.

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

LU C,YUAN L R,GAO S Y, et al. Experimental study on the relationship between tide strength and salt intrusion length [J]. Advances in water science, 2013, 24(2): 251-257.

[11]WANG B, ZHU J, 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): 1 901-1 918.

[12]陈文龙, 邹华志, 董延军. 磨刀门水道咸潮上溯动力特性分析 [J]. 水科学进展, 2014, 25(5): 713-723.

CHEN W L, ZOU H Z, DONG Y J. Hydrodynamic of saltwater intrusion in the Modaomen waterway [J]. Advances in water science, 2014, 25(5): 713-723.

[13]ZHANG Y J, YE F, STANEV E V, et al. Seamless cross-scale modeling with SCHISM [J]. Ocean Modelling, 2016, 102:64-81.

[14]ZHANG Y, BAPTISTA A M. SELFE: a semi-implicit Eulerian-Lagrangian finite-element model for cross-scale ocean circulation [J]. Ocean modelling, 2008, 21(3): 71-96.

[15]ZHANG Y J, ATELJEVICH E, YU H C, et al. A new vertical coordinate system for a 3D unstructured-grid model [J]. Ocean Modelling, 2015, 85:16-31.

[16]ZHANG Y J, STANEV E V, GRASHORN S. Unstructured-grid model for the North Sea and Baltic Sea: validation against observations [J]. Ocean Modelling, 2016, 97: 91-108.

[17]MACCREADY P, GEYER W R. Estuarine salt flux through an isohaline surface [J]. Journal of Geophysical Research Atmospheres, 2001, 106(C6): 11 629-11 638.