强潮河口盐水入侵特征数值模拟研究

路川藤，钱明霞*，夏威夷，丁伟，丁佩

（1.南京水利科学研究院，江苏南京 210029；2.水文水资源与水利工程科学国家重点实验室，江苏南京 210098）

0 引言

入海河口在径流和潮流的共同作用下，形成了特有的物理现象—盐淡水混合。盐淡水混合具有时空分布不均的特点，垂向上，水体可按盐淡水的混合程度分为高度分层型、缓混合型和强混合型[1]；平面上，水体受潮汐强度及海床地形等因素的影响，不同区域盐水浓度不同[2]，同时，不同季节由于径流存在差别，盐水上溯距离也不同[3]。由于影响因素众多，因此盐淡水混合是极其复杂的科学问题。

目前，国内外关于河口盐淡水混合的研究成果主要聚焦在机理和盐水入侵及其工程应用上。在盐淡水混合机理方面，径潮流的此消彼长作用是影响盐淡水平面分布、垂向混合的主要因素[4-5]，此外地形、工程及其气象条件等因素对盐淡水混合也有较大影响[6]；盐水入侵研究表明，除外界因素外，盐淡水混合后导致的水体密度的变化也会对水流流态产生影响[7]，不同河口的影响程度不同；盐水入侵工程应用目前关注较多的是河口水源地饮用水安全问题，例如针对长江口青草沙水库，众多学者在分析研究海平面上升、南水北调、极端天气、人类活动、河势变化等对青草沙水库取水产生影响的基础上，提出了利用三峡水库补淡压咸、北支建闸等工程措施[8-9]，并建立了长江口盐水入侵预报系统[10]，为水源地取水安全提供预警、预报。

本文在前人研究成果的基础上，分析椒江河口盐水入侵特征。椒江位于钱塘江南侧约150 km 处（位置见图1），是我国典型的强潮河口，河口口门处海门站多年平均潮差为4.0 m，最大潮差可达6.3 m[11]。通过本文研究，一方面可以明晰强潮河口的盐水入侵特征，另一方面可以为相关科学研究以及工程建设提供技术支撑和依据。

图1 模型范围与网格示意图Fig.1 Schematic diagram of model domain and grid

1 盐水入侵数学模型构建

1.1 模型范围

本文重点研究椒江河口径流与潮流相互作用下的盐水上溯规律，因此首先建立平面二维盐水入侵数学模型。数值模拟软件采用南京水利科学研究院自主研发的CJK3D-WEM 数值模拟系统[12]，该软件适用于江河湖泊、河口海岸等涉水工程中的水动力、泥沙、水质、温排、溢油的模拟预测。笛卡尔坐标系下的水流以及盐水运动二维方程表达如下：

式中：H为水深，单位：m；z为水位，单位：m；u、v为速度分量，单位：m/s；t为时间，单位：s；f为科氏力，单位：s-1；g为重力加速度，单位：m/s2；Nx、Ny为x、y向水流紊动粘性系数，单位：m2/s；ρ0为水体密度，单位：kg/m3；τsx、τsy为x、y向的风切应力，τbx、τby为x、y向的底切应力，单位：kg/（m·s2）；ϕ为盐度，单位：‰；Kx、Ky为盐水扩散系数，单位：m2/s。

采用三角形网格有限体积法求解式（1），对流项数值通量采用Roe 格式的近似Riemann 解，紊动项采用单元交界面的平均值来计算通过该界面紊动粘性项的数值通量，底坡项采用斜底方法处理，方程的具体求解过程见文献[13]。

数学模型构建以椒灵江干流和台州湾为核心水域，北侧边界至宁波象山县，南侧边界至乐清湾附近，东侧边界至椒江口门以东85 km 处的开阔水域，上游边界包含椒灵江一级支流永安溪和始丰溪（见图1）。模型上游边界采用流量控制，外海边界采用潮位控制。模型网格总数为50 572 个，椒灵江水域网格加密，最小网格边长约为23.5 m，最大网格边长为3 508 m，椒灵江干流河床地形采用2016 年实测地形，台州湾及外海水深数据采用最新海图水深数据，统一到85 高程，坐标系统采用2000 国家大地坐标系。数学模型计算时间步长取1 s；科氏力f=2ωsinφ，其中ω为地球自转角速度，φ为地理纬度，取28.7°；糙率n采用附加糙率方式进行处理，计算公式为：

式中：H为水深。水流紊动粘性系数取0.01HU*（U*为摩阻流速）；动边界水深取0.02 m；盐水扩散系数取30HU*。

1.2 模型验证

数学模型验证包括水动力和盐度两部分。采用2016 年椒灵江河道同步水文测验资料，包括3 个潮位站和5条垂线（见图2），验证时间为2016年12月8日00 时（北京时，下同）—16 日00 时，包括大潮、中潮、小潮的连续过程。图3—4为潮位以及潮流流速流向验证。从图中可以看出，数学模型的潮位计算值与实测值吻合较好，西门、西岙、海门站高低潮位偏差基本在0.1 m 以内，大潮、中潮、小潮的潮流流速及流向与实测值基本吻合，流速偏差基本在10%以内，流向偏差小于15°，数学模型水动力验证总体满足相关规程要求。图5 为各垂线盐度验证，由于实测资料盐度采样时间跨度（3 h）较大，因此盐度的数学模型计算值无法与实际盐度逐时对比，但总体来看，两者大小及过程均吻合良好，盐度自下游向上游逐渐减小，5#点最高盐度值在20‰左右，3#点最高约为12‰，1#点为2.5‰左右。从水动力与盐度的验证结果来看，本文建立的数学模型能够较好地模拟椒江盐水运动特征。

图2 潮位站及潮流垂线位置示意图Fig.2 Schematic diagram of tide level station and tide vertical line position

图3 水位过程验证Fig.3 Tide level verification

图4 潮流过程验证Fig.4 Tidal current verification

图5 盐水入侵过程验证Fig.5 Salinity intrusion verification

2 椒江河口盐水入侵特征研究

2.1 计算条件

椒江为典型的山溪性河流，流量年内分配不均，月度差异明显，根据保证率为50%的流域来水条件分析（见图6），11月平均流量最小，约为23 m3/s，6 月最大，约为684 m3/s，全年平均流量约为160 m3/s，洪枯季流量差异可达30 倍以上。为反映不同径流下的盐水运动特征，上游径流分别取10 m3/s、100 m3/s、200 m3/s、500 m3/s、1 000 m3/s和3 000 m3/s，外海潮位采用1.2 节验证大潮、中潮、小潮潮位过程的结果，海门站大潮平均潮差约为5.2 m，小潮平均潮差约为2.8 m，根据2020 年海门站潮差累积频率统计进行计算，得到大潮累积频率约为5%，小潮约为87%（见图6）。

图6 2001年椒江径流月度分配（a）与2020年海门站潮差累积频率统计（b）Fig.6 Monthly distribution of Jiaojiang River runoff in 2001（a）and statistics of tidal range cumulative frequency of Haimen station in 2020（b）

2.2 盐水分布特征

盐淡水混合是径潮流共同作用的结果。图7与表1 给出了不同径流、不同潮差条件下盐水的最大上溯距离分布（指距离头门岛的距离）。从中可以看到，大潮时盐水上溯距离明显大于小潮，随着径流的增大，盐水逐渐向下游运动，当径流为10 m3/s时，0.1‰盐水可以上溯至始丰溪与永安溪的交汇处三江村附近，5‰盐水可以上溯至八仙岩附近；当径流为200 m3/s 时，0.1‰盐水仅能上溯至西岙附近，距离口外头门岛（见图2中C1处）约60 km，5‰盐水上溯至海门上游13 km 处；当径流为500 m3/s 时，0.1‰盐水仅能上溯至石仙妇上游约6 km 处，距离口外头门岛约40 km，5‰盐水上溯至海门上游4 km处；当径流为3 000 m3/s时，海门以上河段的盐水浓度均小于0.1‰。

表1 不同条件下的盐水最大上溯距离与盐水锋面比降Tab.1 Maximum upstream tracking distance of Brine and the gradient of salt water front under different conditions

图7 不同径流条件下的盐水最大上溯距离Fig.7 Maximum upstream tracking distance of brine under different runoff conditions

在椒江干流沿程布置间距为5～5.5 km 的17个采样点（位置见图2）用于分析不同潮差、不同径流条件下的椒江盐水锋面分布，结果见图8。由图可知，无论是小潮还是大潮，低流量条件下由于盐水上溯距离远，盐水锋面坡度较缓，当径流量为10 m3/s时，小潮的盐水浓度比降约为0.52‰/km，大潮约为0.51‰/km；随着径流量的增大，盐水浓度比降逐渐增大（见表1），当径流量为500 m3/s，小潮的盐水浓度比降约为0.93‰/km，大潮约为0.73‰/km；当径流量为1 000 m3/s，小潮的盐水浓度比降约为0.94‰/km，大潮约为0.76‰/km。综上可见，随着径流量的增大，盐水达到最大上溯距离后，盐水锋面比降变化梯度呈减小趋势。

图8 不同径流潮差条件下的盐水锋面分布Fig.8 Salt water front distribution under different runoff and tidal ranges

2.3 径流对盐水上溯距离变化梯度影响分析

在特定的径流条件下，盐水在涨落潮的作用下最大上溯位置不断发生变化，涨憩时盐水上溯距离最远，落憩时最近。图9 统计了椒江河口盐水在大潮和小潮不同潮型作用下的盐水上溯距离变化梯度与径流变化的关系，即：

图9 盐水上溯距离变化与径流的关系Fig.9 Relationship between saltwater upstream tracking distance change and runoff

式中：ε 为盐水上溯距离随径流的变化梯度，单位：m/（m3·s）；∆Q为不同径流Q1、Q2的差值，单位：m3/s；∆d为径流Q1、Q2下盐水上溯距离d1、d2之间的差值，单位：m。

由图9可以看到，无论在涨憩还是落憩时刻，随着径流量的增大，盐水上溯距离呈减速降低趋势，涨憩时刻盐水上溯距离变化梯度略大于落憩，小潮盐水上溯距离变化梯度大于大潮。在10～100 m3/s的流量区间，涨憩时刻大潮盐水上溯距离变化梯度为222 m/（m3·s），小潮为238 m/（m3·s），落憩时刻两者数值分别为207 m/（m3·s）和233 m/（m3·s），而在200～500 m3/s流量区间，盐水上溯距离变化梯度急剧降低至15～35 m/（m3·s），这说明在径流抑制盐水上溯的过程中，盐水锋面距离河口越远，其受径流的影响越大，反之，受径流影响相对较弱。

2.4 潮差对盐水上溯变化区间影响分析

根据数学模型计算结果，统计了在特定径流条件下的小潮、中潮和大潮（对应潮型25 h 内的最大潮差分别为3.2 m、5.1 m、6.2 m）盐水上溯距离变化区间（涨憩盐水上溯距离-落憩盐水上溯距离），结果见图10。由图可知，随着潮差的增大，盐水上溯变化区间总体呈线性增长趋势，随着径流量的增大，盐水上溯变化区间逐渐减小。当径流为10 m3/s时，小潮盐水上溯变化区间约为12 km，大潮约为21 km，而径流为500 m3/s 时，大小潮的盐水上溯变化区间明显下降，小潮约为8 km，大潮约为14 km，说明盐水锋面距离河口越远，潮差对盐水上溯变化区间影响越大，反之，盐水锋面距离河口越近，潮差对盐水上溯变化区间影响越小。

图10 潮差对盐水上溯区间影响Fig.10 Influence of tidal range on salt water upstream tracking region

2.5 椒江河口盐水上溯距离函数关系式

自然条件下，盐水上溯距离主要与径流和潮差有关。根据表1中径流与盐水上溯距离的关系分析可知，在潮差不变的条件下，盐水上溯距离与径流的自然对数（lnQ）具有一定线性关系，因此用lnQ代表径流因子，同时在径流不变的条件下，盐水上溯距离与潮差的线性关系较好（见图10）。由以上分析可以构建椒江河口盐水最大上溯距离与径流、潮差的函数关系式：

式中：d为盐水上溯距离（即距离头门岛的距离，单位：km）；Q为径流，单位：m3/s；∆h为潮差，单位：m；a、b、c为系数，由于盐水上溯距离与径流呈反比，与潮差呈正比，因此a＜0，b＞0。

根据数学模型计算不同径流、不同潮差下涨憩和落憩各18 组盐水上溯的最大距离，利用式（4）进行拟合，求解a、b、c。经拟合计算，涨憩时a= -9.95，b= 2.18，c= 89.67；落憩时a= -8.08，b= 0.34，c= 75.96。图11为盐水上溯距离函数关系式预测值与数学模拟值的对比结果，经统计，涨憩时刻两者偏差为2%～3%，落憩时刻为5%～6%。涨憩和落憩时刻数学模拟值与预测值的相关性高达0.99，说明本文建立的椒江河口盐水上溯距离函数关系式预测精度相对较高，可为椒江盐水相关研究提供借鉴。

图11 盐水上溯距离函数关系式预测值与数学模拟值对比Fig.11 Comparison between predicted value and digital analog value of salt water upstream tracking distance function

3 结语

以强潮河口——椒江河口为研究对象，通过数学模型研究了不同径流、潮差条件下的盐水入侵特征。主要结论如下：

①椒江河口是典型的强潮河口，当径流为10 m3/s 时，0.1‰盐水可以上溯至始丰溪与永安溪的交汇处三江村附近，5‰盐水可以上溯至八仙岩附近，当径流为500 m3/s 时，0.1‰盐水仅能上溯至石仙妇上游约6 km 处，5‰盐水上溯至海门上游4 km 处，当径流为3 000 m3/s时，海门以上椒江河段的盐水浓度均小于0.1‰。

②在径流变化条件下，盐水上溯距离变化梯度涨憩大于落憩，小潮大于大潮，随着径流量的增大，盐水上溯距离变化梯度逐渐增大；在潮差变化条件下，盐水上溯距离变化区间呈线性增长趋势，随着径流量的增大，同等潮差下的盐水上溯距离变化区间逐渐减小。

③通过数学模拟结果拟合出了椒江河口盐水最大上溯距离与径流、潮差的函数关系式，涨憩时刻关系式预测值偏差为2%～3%，落憩时刻为5%～6%，预测值精度较高，可为椒江盐水相关研究提供借鉴。