盐水入侵是河口的一种普遍现象，盐水入侵将改变水体的理化特性，对河口水流结构、泥沙输运有重要影响，同时对淡水资源利用来说是一种海洋灾害^[1]。因此研究盐水入侵既具有重要的科学价值，又具有实际应用价值。

长江口年径流量大，潮汐中等，其盐淡水混合特征属于缓混合型。河口盐水入侵主要受径流量和潮汐的影响^[2]，还受地形、风应力、科氏力以及口外陆架环流等的作用^[3]。长三角地区经济发达，人口稠密，需要大量淡水资源，长江口有4个主要水库：陈行、宝钢、青草沙以及东风西沙(见图 1)。当水库盐度超过0.25 g/L时(约0.45 psu)，水源地应停止取水。2014年中国海平面公报统计显示：1980—2014年中国沿海海平面平均上升速率为3.0 mm/a^[4]。海平面上升将导致更多高密度盐水入侵河口地区，河口淡水资源将受到更严重的威胁。国外学者就海平面上升对河口区域的盐水入侵进行了相关研究：海平面上升1 m可以使孟加拉国的Gorai河盐度增加约1.5 psu^[5]，使美国的James河枯季取水口的平均盐度增加5 psu^[6]。随着海平面上升，美国Chesapeake湾的平均盐度、盐水入侵距离以及盐度分层都将增加^[7]。许多学者对长江口盐水入侵方面做了深入的研究：朱建荣等研究了陈行^[8]、青草沙^[9]和东风西沙^[10]水库的盐水入侵和人类活动(三峡工程^[11]、南汇边滩围垦工程^[12]以及北支新村沙围垦^[13])对长江口盐水入侵的影响；张二凤等^[14]采用实测资料研究北支盐水入侵的规律及影响因素，研究表明径流量小、潮差大以及强偏北风三者的叠加是导致强盐水入侵的主要原因；戴志军等^[15]研究了2006年特枯水文年汛期情景对淡水资源的影响，研究结果表明：北支倒灌南支严重，导致陈行、宝钢水库不存在淡水资源。结合海平面上升对长江口盐水入侵方面的研究比较少：杨桂山等^[16]通过相关性分析法分析了海平面上升对吴淞和高桥盐度的影响，分析结果表明海平面上升0.8 m情形下，只有当径流量达到15 000 m³/s时盐度才不受海平面上升的影响。裘诚^[3]采用通量机制分解方法分析了海平面上升对长江口盐水入侵的影响，研究表明海平面上升导致垂向环流加强，增强了垂向分层结构。

然而海平面上升对长江口盐水楔的影响，至今未见报道。因此本文基于MIKE3建立了长江口三维水动力及盐度输运数学模型，分析长江口枯季盐水入侵对海平面上升的响应特征。

1 模型建立与验证

丹麦水力学研究所(DHI)研发的MIKE3模型，主要模拟海洋、河口、河流、湖泊及海岸的水流、波浪、泥沙及环境变化，为海洋管理和规划以及工程应用提供了完备、快捷、有效的设计环境。MIKE3 FM属三维潮流模型，根据静水压力假设、Boussinesq假设、浅水条件和适定边界条件，采用控制体积法求解Navier-Stokes方程概化的浅水方程^[17]。

1.1 模型区域和网格

为了更好地模拟冲淡水的影响，模型模拟范围东至长江口口外200多千米，北到吕泗，南到象山以南，西到长江江阴以及钱塘江仓前，包括整个长江口和杭州湾区域(图 2)。长江口三级分汊、四口入海的格局以及舟山众多岛屿导致岸线特别曲折，三角形网格更贴合岸线，因此模型采用SMS构造三角形网格，模型共有18 006个节点，33 656个单元。空间步长在130~29 880 m的范围内。

1.2 边界条件和参数设置

模型中的底部切应力由公式${\tau _{\rm{b}}} = \rho g\left| {{u_{\rm{b}}}} \right|{u_{\rm{b}}}/\left( {\frac{1}{\kappa }\ln \left( {\frac{{\mathit{\Delta }{z_{\rm{b}}}}}{{m{k_{\rm{s}}}}}} \right)} \right) $计算，其中：κ为冯卡门常数，取0.4，m为系数取1/30，k_s为糙率高度，其值为0.000 27~0.001 6 m，u_b为距离海床Δz_b处的流速，垂向采用σ坐标，均匀分为5层。水平扩散和垂直扩散均采用scaled eddy viscosity公式，水平扩散比例因子为1，垂直扩散比例因子为0.01。紊流采用标准k-ε模型计算，垂直涡黏性系数υ_t=c_μk²/ε，经验系数c_μ取0.09，k为紊动能，ε为紊动能耗散率。模型采用干湿动边界处理技术，模型中干点临界水深取0.005 m，湿点临界水深取0.05 m。模型的计算时间步长在0.1~30.0 s的范围内自行调节，从而保证模型计算过程中始终满足克朗数(CFL)小于1。模型的水平涡黏性系数采用Samagorinsky亚网格尺度模型来计算，其中Samagorinsky系数取为常数0.28。初始流速和潮位分别为1 m/s和0.5 m。盐度初始采用计算稳定后(冷启动，6个月后)的场值，江阴和仓前的盐度取为0，外海边界盐度根据实测资料插值给定。

模型验证时江阴边界采用大通站实测流量过程控制，仓前取多年平均流量1 000 m³/s^[18]，风应力采用美国国家海洋和大气管理局(NOAA)提供的海面10 m以上处测量间隔为6 h的风场资料来计算。

1.3 模型验证

采用长江口实测潮位、流速、流向以及盐度资料对模型进行了验证，潮位实测资料时间为2005年8月14—24日^[19]；流速流向实测资料时间为2004年5月5—6日^[19]；盐度实测资料时间为2003年2月17日11：00(大潮)—18日以及2月23日14：00—24日13：00(小潮)^[20]。实测站点见图 1，潮位、流速流向以及盐度验证过程分别见图 3~5。验证结果表明，潮位计算值与实测值基本吻合，潮流计算值在大小及相位上与实测值大体一致，模拟的盐度过程基本能体现涨落潮规律。因此本模型可以用于计算分析海平面上升对长江口盐水入侵的影响。

2 径流量和海平面上升值的选取

预测分析海平面上升对长江口盐水入侵影响时，为了模拟长江口盐水入侵最严重情形，绘制1950—2011年大通1月流量累计频率曲线，径流量的丰、平、少年际变化是根据流量点绘累积频率曲线来判断，频率小于25%为丰水年，大于75%为少水年，25%~75%之间为平水年，因此江阴采用大通枯季1月75%频率的流量，即8 900 m³/s，仓前边界采用钱塘江多年平均枯季流量600 m³/s。

海平面上升是一种缓发性灾害，其长期的累积效应会淹没滨海低地等，对沿海地区的经济发展构成威胁。海平面上升速率主要根据验潮站和卫星高度计资料获得，因此海平面上升速率随着研究区域及时间范围不一样相差较大。2007年第4次IPCC报告指出^[21]，1961—2003年期间的验潮站资料分析显示全球海平面上升平均速率为1.3~2.3 mm/a，而1993—2003年平均上升速率为3.1 mm/a。2013年第五次IPCC报告指出^[22]，1901—2010年全球海平面上升速率为1.7 mm/a。2014年中国海平面公报显示：1980—2014年中国沿海海平面上升速率为3.0 mm/a。研究海平面上升对河口盐水入侵的影响，需要预测未来实际的河势情况，要给出相对海平面上升，即在绝对海平面上升的基础上，还需要考虑地面沉降。刘毅和龚士良^[23]在上海地面沉降长期预测研究的实际经验基础上，预测长江口未来地面沉降速率在1.82~6.73 mm/a。为了研究海平面上升最不利的影响，按海平面上升速率3.1 mm/a和地面沉降速度6.73 mm/a计算未来100年后相对海平面可能上升值，即约1 m。因此，本文计算原海平面和海平面上升1 m情形下的长江口盐度输运，以此分析海平面上升对长江口盐水入侵的影响。为了取得稳定值，模型采用热启动，模型模拟时间为2004年10月1日到2005年1月31日，平均盐度取2005年1月份的计算结果进行分析，选取小潮涨急时刻的盐度垂直结构进行分析。

3 海平面上升引起的盐度平面变化

河口地区的盐水随涨潮流向上游运动，随落潮流向口外运动，往往涨憩时刻盐水向上游入侵最远，落憩时刻向口外海岸退得最远。图 6为原海平面及上升1 m情形下涨憩时刻盐度的平面分布。原海平面情形下，南北支交界处，北支盐水倒灌南支，使得崇头出现小部分高盐水水体(0.45~2.00 psu)，南支的月平均盐度大部分都小于0.45 psu，未超过国家饮用水的标准，适合取水。北港的月平均盐度相对较高，为0.45~2.00 psu，横沙通道是北港和北槽水体交换的重要通道，北港盐水入侵除来自自身通道外还来自横沙通道，因此北港南部盐度比北部高。

海平面上升1 m情形下，北支倒灌进入南支的水体盐度增大，而且影响范围也增大，南支水域的平均盐度均超过0.45 psu，南港水域的平均盐度为2.00~20.00 psu，北港水域的盐度为0.45~10.00 psu。北支中下段的盐度减小，是由于海平面上升使得北支分流比增大^[24]，北支下段常年被高盐水包围，涨潮动力的增强对高盐水包围的北支下段影响较小，反而北支落潮分流比的增加使得北支中下段盐度减小。九段沙、横沙东滩以及崇明东滩地区滩地边缘高程小于1 m区域被淹没，使得滩地被盐水覆盖的面积增大。

4 盐水楔对海平面上升的响应

当径流注入中等或弱潮河口时，由于海水密度较淡水密度大，淡水将从表层泄入海中，海水从底部侵入，呈“楔状”，故称为盐水楔。选取枯季小潮涨急时刻盐度垂直结构分析海平面上升对长江口盐水楔的影响。涨急时刻以图 6(a)中五角星位置处为参考点。为了进一步了解海平面上升引起的湍流混合和层化现象的变化规律，可通过“梯度Richardson数”来衡量：${R_i} =-\frac{g}{\rho }\frac{{\partial \rho }}{{\partial z}}{\left( {\frac{{\partial u}}{{\partial z}}} \right)^{-2}} $，其中$\frac{{\partial \rho }}{{\partial z}} $和${\frac{{\partial u}}{{\partial z}}} $分别为某一水层的密度差(kg/m³)和流速差(m/s)，ρ为垂向平均密度(kg/m³)，∂z为某一水层的水深(m)。数学模型垂直分为5层，垂直间隔比较大，计算的Richardson数不能反映小于1 m水层中的混合程度，本文根据文献[25]采用比理论Richardson数0.25更大的临界值1来判别海平面上升对水体的层化和混合影响，当Richardson数小于1时，水体是混合的；当Richardson数大于1时，水体出现层化现象，且Richardson数越大，层化现象越明显。

在北支和南支深槽布置两条纵剖面D1和D2(位置见图 1)，图 7为原海平面及上升1 m情形下北支下段河道纵剖面(D1)上的盐度垂向分布，由于北支的径流分流比非常小(小于1%)，导致盐水入侵严重，盐度比较高，还将倒灌南支。图 8为原海平面及上升1 m情形下北支下段AA1点(位置见图 7(a))的Richardson数垂向分布。水体中是否会产生层化的现象，取决于河口环流和潮汐应变的能量和潮汐搅动能量是否相互平衡。原海平面情形下，北支有部分径流注入，出现盐水楔现象，表层Richardson数为2左右(大于1)，中间层Richardson数达到60，而底层小于1，说明中上层有层化现象，且中层远大于表层，底层则是混合较好；海平面上升1 m情形下，北支下段等盐线更稀疏，盐水楔变缓并减弱，层化现象上移，表层Richardson数达150。海平面上升一定高度使得北支盐水楔现象增强，主要原因是因为海平面上升增加了北支上段的落潮流，使得表层盐度减小，从而加剧北支下段的盐水楔。

图 9为原海平面及上升1 m情形下南支、南港和南槽河道纵剖面(D2)上的盐度垂向分布，南支和南港不与外海直接相连，且主要受淡水的控制，盐度均小于0.45 psu，盐度表底层差异不明显。

图 10为原海平面及上升1 m情形下AA2点的Richardson数垂向分布。原海平面情形下，在IP1~IP18段，垂向分层不明显，无明显盐水楔发育，盐水楔出现在IP19~IP22段，即南槽的中下段，而在IP22~IP23段，由于径流作用较弱，表底层水体盐度无差异，也无明显盐水楔发育，AA2点的Richardson数说明层化现象上层＞底层＞中层；海平面上升1 m情形下，盐水楔整体向上游移动较小距离，等盐线更密集，AA2点的Richardson数底部减小约2~4，顶部增大16。

5 结语

本文利用MIKE3软件建立了长江口三维水动力及盐度输运数学模型，分析盐水入侵最不利的影响，基于1950—2011年大通1月流量累计频率曲线，采用大通枯季1月75%频率的流量8 900 m³/s，计算了海平面上升对长江口枯季盐水入侵的影响，得到了如下主要结论：

(1) 海平面上升后，北支上段、南支、南北港以及南北槽的盐度都明显上升，南北支水体交换增强，更多南支低盐水进入北支中下段，使得北支中下段盐度明显减小；北支倒灌进入南支的高盐水更多使得南支盐度增大。

(2) 海平面上升1 m后，北支下段盐度层化现象减弱并由中间层移向表层，表层Richardson数达150，这是由于潮汐增强，搅动作用增强，层化现象减弱，但是北支径流增大，表层盐度减小，会使得表层层化增强。

(3) 海平面上升1 m后，南槽盐水楔上移较小距离，底部Richardson数减小2~4，顶部增大16，但垂向变化趋势不明显。