钱塘江河口盐水入侵影响因素及数值预报

潘存鸿，张舒羽，史英标，鲁海燕

(1．浙江省水利河口研究院，浙江 杭州 310020;2．浙江省海洋规划设计研究院，浙江杭州 310020)

盐水入侵是潮汐河口存在的普遍问题，研究盐水入侵对河口水动力、泥沙运动、河床冲淤、水资源利用、河口环境生态等都具有重要的学术意义和实际应用价值，一直是国内外研究的热点［1-6］。水流氯度数值模拟是研究盐水入侵规律，流域建库、河口治理、跨流域引水等人类活动对河口咸水入侵影响，以及预测、预报江水氯度的重要研究手段［1，7-10］。

杭州湾是世界上著名的强潮河口湾，湾顶澉浦站多年平均潮差5.57 m，历史最大潮差9.0 m。乍浦以上河床溯源抬高，水深减小，潮波变形显著，在澉浦上游形成世界闻名的钱塘江涌潮。钱塘江河口涨潮量大，而径流相对较小，富春江电站断面多年平均径流量952 m3/s，且年际、年内分布不均。因此，钱塘江河口盐水入侵长度和强度远大于其他潮汐河口。在自然径流条件下，盐水可上溯到离杭州湾口207 km处的闻堰以上。目前，杭州市区自来水厂的水量85%取自钱塘江河口闸口段(图1)，枯水大潮期，咸潮倒灌，不同程度地影响了杭州市自来水厂取水口的水质，如在1978年、1989—1995年、2003年等年份出现了比较严重的咸水入侵，导致自来水厂取水口连续数天不能取水，严重影响了杭州的工农业生产和市民的正常生活。早在1981年韩曾萃等［11］应用一维水流氯度数值模型预测了杭州市钱塘江取水口处水体氯度，2012年韩曾萃等［12］又建立了钱塘江河口一维水流氯度动床数值模型，结合实测资料统计分析，用于指导富春江水库每天的拒咸下泄径流量调度，以保证杭州市供水不受盐水入侵的影响。

图1 钱塘江河口形势

本文在分析影响钱塘江河口盐水入侵的主要因素的基础上，建立了考虑涌潮作用的钱塘江二维氯度预报数值模型，以预报短期(1个月以内)、中期(1～3个月)杭州市公共水厂钱塘江取水口处的氯度及超标时间。

1 盐水入侵影响因素分析

钱塘江河口山潮水比值仅为0.02［8］，属强混合型河口，盐水入侵严重。影响钱塘江河口盐水入侵的因素众多，除一般潮汐河口所共有的径流、潮汐、风速风向等影响因素外，江道地形、取水流量和杭州湾氯度大小等对钱塘江盐水入侵也有较大影响。鉴于盐度S与氯度C存在很好的线性关系(即S=1.80655C)，因此本文采用氯度代表盐度进行分析。钱塘江河口闸口、七堡、仓前、盐官、澉浦等水文站(位置见图1)每日在日潮低平潮和高平潮时取样测量氯度，分别代表潮周期内氯度的最小值和最大值。本文基于这些实测资料，结合数学模型成果，重点分析径流、潮汐(包括涌潮)和江道地形对氯度变化的影响。

1．1 径流对盐水入侵的影响

径流是影响钱塘江河口盐水入侵的最重要因素之一，径流变化及其对盐水入侵影响的特点如下:

a．径流年际变化对盐水入侵的影响。钱塘江河口多年平均流量为952 m3/s，最大年平均流量为1710 m3/s(1954年)，最小年平均流量为411 m3/s(1979年)，流量年际间变幅大。年际间有连续丰水年和连续枯水年的周期变化［13］，20世纪60年代和80年代(除1983年)及21世纪前10年，径流偏枯，20世纪50年代和70年代及90年代则偏丰。由1994—2013年年径流量与七堡站年平均氯度的关系(图2)可以看出，由于年内径流分配差异较大，且存在江道地形冲淤引起的潮汐变化，导致年径流量与年平均氯度相关性不好，说明不宜采用年尺度来分析径流与氯度的关系。

图2 径流量与氯度关系曲线

b．径流年内变化对盐水入侵的影响。钱塘江河口径流年内变化不均，存在明显的丰水期和枯水期变化，丰水期4个月(3—6月或4—7月)径流量占全年的70%左右［14］。径流的季节性变化特点导致盐水入侵的年内变化，在丰水期盐水入侵距离短，沿程江水氯度低;反之，在枯水期(特别是秋季大潮枯水期)盐水入侵严重，沿程江水氯度高。由1994—2013年丰水期 (4—7月)4个月径流量与相应时间的七堡站平均氯度的关系(图2)可知，氯度与径流量基本上成反向变化趋势，径流量越大，氯度越低。

c．受上游建库影响显著。1949年以来，钱塘江流域兴建了小型以上水库 2000余座，其中1960年建成的新安江水库(总库容为216亿m3)对径流的影响尤为明显。水库建成后丰水期径流量减少，枯水期径流量增大，径流量趋于均匀。另外，建库还削减了洪峰流量，洪峰越大削减越多;同时，还减少了大洪水的出现次数。对杭州河段盐水入侵而言，上游建库引起的径流调节作用主要表现在枯水期流量的增大，相应地咸水入侵距离减小了20 km左右［8］。

图3为潮汐条件不变(七堡潮差2.8 m)，不同径流条件下钱塘江河口咸水(以250 mg/L为标准)的入侵距离，可见，随着径流增加，咸水入侵上溯距离减小，且递增相同流量的情况下，本底流量越小，上溯距离减小越明显。

3 不同径流条件下钱塘江河口咸水入侵距离

1．2 潮汐及涌潮对盐水入侵的影响

潮汐是影响钱塘江河口盐水入侵的另一个重要因素。钱塘江河口属正规半日潮，一日两潮，并存在半个月的大小潮变化周期。图4为盐官站2012年10月9—24日半个月实测氯度和潮位变化，其中氯度仅在日潮取样测量2次，低平潮和高平潮时取样，分别代表潮周期内氯度的最小值和最大值。由图4可知，氯度变化基本与潮汐变化相同，大潮期间无论是最大氯度还是最小氯度均较高，而小潮期间最大氯度和最小氯度均较低，氯度变化也存在明显的半月潮汛变化周期。

图4 盐官站实测潮位和氯度

图5为盐官站1个潮周期内潮位、流速和氯度同步变化过程线，可明显看出，在高潮位后涨憩附近氯度达到最高，在低潮位落憩时刻(涌潮前)氯度达到最低。在相位上，氯度变化略滞后于潮汐变化。

钱塘江以涌潮闻名于世。大潮时，钱塘江涌潮在澉浦上游开始形成，向上游传播过程中逐渐增大，至大缺口—盐官河段达到最大(潮头高度最大可达3 m以上)，以后逐渐减弱，最远可上溯至闻家堰以上，涌潮河段超过100 km。涌潮对盐水入侵的影响主要表现在［15］:①涌潮促使盐水入侵加剧。对固定位置而言，一般潮差越大，涌潮越强，盐水入侵加剧。②涌潮过后纵向上氯度梯度增大。涌潮前后，潮位、流速发生突变，从而使涌潮前后氯度也发生剧变，如图5所示。

图5 盐官站潮位、流速和氯度过程线

1．3 江道地形对盐水入侵的影响

钱塘江河口潮强流急，涌潮汹涌，径流变幅大，泥沙易冲易淤，在径流和潮流作用下，河床存在冲淤幅度大、速度快等特点［14］。因钱塘江河口河床宽浅，河床冲淤变化反过来又影响潮汐大小。一般钱塘江河口大缺口以上具有洪冲潮淤的特点，年际丰水年江道容积大，枯水年江道容积小;年内每年汛后8月江道容积最大，以后逐渐淤积，至次年2月江道容积最小。表1为江道容积较大年份(丰水年1995年)和较小年份(枯水年2004年)的中潮位下江道容积、七堡月平均潮差和七堡月平均氯度的比较。由表1可知，1995年8月江道容积、七堡月平均潮差、七堡月平均氯度分别是2004年同期的3.4倍、2.8倍和2.1倍;1995年11月江道容积、七堡月平均潮差、七堡月平均氯度分别是2004年同期的2.6倍、7.2倍和19.4倍。因此，江道容积越大，低潮位下降越明显，且阻力减小越多，从而潮差和进潮量越大，导致盐水入侵越严重。

表1 江道地形对盐水入侵的影响

钱塘江河口大缺口以上汛后江道地形主要取决于汛期的径流量，径流量越大，汛后的江道容积越大［14］。从这一角度看，江道地形对盐水入侵的影响也可以认为是径流对盐水入侵的间接影响，是通过河床地形冲淤反馈变化引起潮汐变化，从而间接影响盐水入侵的，在时间上存在一定的滞后性。

澉浦以下的杭州湾因水深大，河床冲淤对潮汐变化影响极小，大缺口以上因水深小，河床冲淤对其上游的潮汐变化影响较大。因此，当选用澉浦以下的潮汐作为影响盐水入侵的因素时，需考虑江道地形的变化;若选用上游河段的潮汐(如七堡站)则可以不考虑江道地形的影响，因为上游河段的潮汐已包含了下游潮汐和江道地形的综合影响。

2 二维氯度预报数值模型的建立和验证

钱塘江河口河床宽浅，潮强流急，涌潮汹涌，水流充分混合，与一般潮汐河口相比，氯度在垂向上分布较为均匀，但在横向上仍存在一定差异，而三维数学模型计算量大，在时间上难以满足实时预报要求。为此，建立考虑涌潮作用的二维氯度预报数值模型来预测钱塘江河口氯度的时空变化。

2．1 模型建立

在涌潮水流数学模型的基础上［16］，增加了氯度对流扩散方程，同时还在动量方程中考虑了氯度密度引起的压力项，以及取排水引起的水量、动量和氯度的变化，来建立考虑涌潮作用的二维氯度预报数值模型。模型控制方程为非恒定二维浅水水流和氯度对流扩散方程的守恒形式:

式中:u、v分别为x、y方向的流速分量;h为水深;g为重力加速度;C为氯度;Cs为排水口的氯度;Sbx、Sby分别为 x、y 方向的底坡源项;Sfx、Sfy分别为 x、y 方向的阻力项:Ex、Ey分别为 x、y方向的扩散系数;ρ为盐水密度;Qs为取(排)水口的水流流量;A为取排水口位置的计算单元面积;β为取(排)水方向与x轴的夹角;b为床底高程;n为Manning系数。

计算域采用三角形网格离散，并采用单元中心格式。设Ωi为第i个三角形单元域，Γi为其边界，对式(1)应用有限体积法离散，经推导化简可得

式中:Ai为三角形单元 Ωi的面积;(cosθ，sinθ)为 Γ外法向单位向量;Δt为时间步长;lj为三角形边长;下标j表示i单元第j边;上标n为时间步。

采用有限体积法求解式(2)，其中水流法向数值通量采用能模拟涌潮间断流动的KFVS格式求解，详见文献［16］。氯度法向数值通量采用基于准确Riemann解的Godunov格式求解，该格式能模拟浓度的强间断问题。设单元i与k的共边(界面)为Γik，引进局部坐标ξ，ξ与Γik垂直，且从i单元指向k单元。略去源项和二阶项后，氯度对流扩散方程对应的Riemann问题可表达为［17］

其中 Ec=hC Fc=hωC

式中:ω为ξ方向的流速，即法向流速;Ecl和Ecr为常数;下标l、r表示求解通量的界面的左侧和右侧。

式(3)有4种可能的解:常数解、稀疏波、激波、滑移线。后3种解为:①稀疏波。稀疏波扇区内氯度C与初始状态C0的关系为C=C0。②激波。根据Rankine-Hugoniot间断条件，可得氯度C=C0。③滑移线。根据Rankine-Hugoniot间断条件，可得氯度

应用上述常数解、稀疏波、激波及滑移线关系式求得Riemann解后，可得界面的法向数值通量。采用与MUSCL类似的方法，构建了空间二阶精度的计算格式［18］，以提高计算精度。

2．2 模型验证

模型验证时计算下边界为澉浦—西三断面，计算上边界分别为富春江水电站和浦阳江湄池，计算区域采用三角形网格剖分。选取2007年和2011年秋季1个月左右的水流和氯度资料作为模型验证资料。计算地形采用相应年份的7月或11月1/5万测图。水流边界条件为:下边界给定实测潮位过程，两个上边界给定流量过程。氯度边界条件为:流入时给定氯度值，流出时采用纯对流方程计算得到。

因氯度计算的初始条件对计算结果影响的时间较长，为此，将验证计算时间提早1～2个月开始，以消除氯度初始条件对计算结果的影响。计算验证了沿程盐官、仓前、七堡、闸口、闻家堰等站的潮位和氯度资料，以及水文测验期间的临时潮位、流速、涌潮等资料。图6和图7分别为2011年9月仓前站的潮位和氯度验证结果，计算值与实测值基本吻合;图8为仓前站的涌潮验证结果，图中涌潮到达后实测潮位数据为每分钟1个值，可见涌潮高度计算值与实测值相当接近。

图6 仓前潮位验证

图7 仓前氯度验证

图8 仓前涌潮验证

模型计算结果表明:涌潮促使盐水入侵加剧;涌潮过后氯度的时间梯度和空间梯度迅速增大，涌潮是形成氯度锋的动力机理。上述计算结果与实测资料分析是一致的。

3 二维氯度预报数值模型的应用

3．1 计算条件

建立的钱塘江河口二维氯度预报数值模型不但可用于研究盐水入侵规律、氯度时空分布规律、人类活动对盐水入侵的影响等，并且还可用于氯度预报。考虑到模型为定床模型，而钱塘江河口冲淤复杂，一般不适用于3个月以上的长期预报，仅适用于短期(1个月以内)、中期(1～3个月)预报，包括实时预报。由于预报时一般只给定上游富春江电站径流量，需预报时刻的潮汐、江道、下边界氯度等条件未知，因此需分析计算条件。

a．上、下边界水流条件:对澉浦站进行调和分析，利用调和常数预报得到逐时的天文潮潮位过程，作为下边界澉浦的水边界条件;浦阳江上边界给定枯水流量;富春江电站给定径流量。

b．江道地形条件:钱塘江河口澉浦至闸口每年有4月、7月和11月3次江道地形测量图，氯度预报一般预报盐水入侵较为严重的秋季(8—11月)枯水大潮期，因此，可采用当年的7月测图作为计算地形。而闸口以上江道地形变幅较小，且对计算结果影响不大，尽量采用较新的地形图即可。

c．下边界澉浦断面氯度:澉浦断面氯度不但受上游钱塘江流域径流丰枯的影响，并且还受下游杭州湾水体氯度的影响，而杭州湾水体氯度则受到长江口前期径流的影响。根据1971—2010年历年8—11月澉浦站氯度统计结果，8—11月每天最大氯度平均值为4952 mg/L，结合2011年水文测验期间澉浦断面氯度实测资料，选取5000 mg/L作为涨潮时氯度边界条件。

d．杭州河段取排水流量较大，除公共水厂总取水流量25.7 m3/s和七格污水处理厂、四堡污水处理厂合计排污流量13.3 m3/s变化不大外，其他两岸取排水流量有较大的不确定性。根据调查分析，模型另增加两岸引水流量50～100 m3/s。

应用模型进行氯度预报时，应先采用预报当年前期的实测潮汐、氯度资料进行检验，检验取排水流量、下边界澉浦断面氯度选用的合理性，必要时可根据前期实测资料和调研情况进行修正，以提高模型的预报精度。

3．2 模型预报结果

因氯度计算的初始条件对计算结果影响的时间很长，为让模型适应实时预报，必需提高模型计算效率。除在计算区域选取、网格布置上考虑计算机时外，还特别采用了 OpenMP并行编程技术［19］，对基于多核处理器的计算程序实现了并行计算，以加快计算速度。经比较，在4核处理器的计算机上运行，计算时间可缩短2/3左右，计算效率大大提高，模拟30 d盐水入侵过程耗机时仅30 min左右(CPU Intel Xeon(R)E31225，主频 3.1 GHz，内存 4 G)。

应用二维氯度预报数值模型对2012年10月钱塘江河口杭州市公共水厂取水口氯度进行了预测。模型采用的江道地形闸口至澉浦选取2012年7月测图，闸口至闻家堰选取2012年2月测图，闻家堰以上选取2011年7月测图。水流边界条件:下边界澉浦断面给定潮汐预报的天文潮潮位过程;上边界富春江电站给定日平均流量。氯度边界条件:下边界澉浦断面流入时给定2011年水文测验期间实测平均氯度 5000 mg/L，流出时由纯对流方程求得。杭州河段取水流量采用110 m3/s、排污流量采用13.3 m3/s。预报结果与后期获得的实测资料比较见表2，因实测资料仅当氯度大于250 mg/L才有记录，表中仅列出了超标的最大氯度和超标时间。由表2可知，无论是最大氯度还是超标时间，预报结果与实测值基本吻合。

表2 南星水厂取水口2012年10月14—20日最大氯度和超标时间预报与实测结果

此外，为适应非专业人士采用本模型，将二维氯度预报数值模型集成到钱塘江河口预报预警系统中，该系统由信息管理、氯度预报、氯度预警等功能模块组成［20］。

4 结论

a．影响钱塘江河口盐水入侵的主要因素有径流、潮汐及涌潮、江道地形等。不同于其他潮汐河口，因钱塘江河口河床冲淤幅度大，水深浅，河床冲淤对潮汐大小影响很大。江道容积越大，低潮位下降越明显，且阻力减小越多，从而潮差和进潮量越大，导致盐水入侵越严重，江道地形通过改变潮汐大小间接影响盐水入侵，其在时间上存在滞后，是钱塘江河口所独有的特性。

b．涌潮促使盐水入侵加剧;涌潮过后氯度梯度在空间上和时间上均明显增大，涌潮前后氯度存在突变，是形成氯度锋的动力机理。

c．在二维涌潮水流数值模型的基础上，应用基于准确Riemann解的有限体积法建立了具有空间二阶精度的二维氯度预报数值模型，可用于钱塘江河口短、中期氯度实时预报，模型具有较大的普适性和实用价值。采用OpenMP并行编程技术大大提高了计算效率，使模型适用于氯度实时预报;同时将氯度预报模型集成到钱塘江河口氯度预报预警系统，以方便非专业人士使用。

［1］沈焕庭，茅志昌，朱建荣．长江河口盐水入侵［M］．北京:海洋出版社，2003．

［2］SAVENIJE HHG．Salinity and tides in alluvial estuaries［M］．Amsterdam:Elsevier Press，2005．

［3］陈水森，方立刚，李宏丽，等．珠江口咸潮入侵分析与经验模型:以磨刀门水道为例［J］．水科学进展，2007，18(5):751-755．(CHEN Shuisen，FANG Ligang，LI Hongli，et al．Saltwater intrusion analysis and experiential model for Pearl River estuary，South China:a case study in Modaomen Watercourse［J］．Advances in Water Science，2007，18(5):751-755．(in Chinese))

［4］DUC N A．Salt intrusion，tides and mixing in multi-channel estuaries［M］．Amsterdam:Taylor ＆ Francis Group，2008．

［5］ROSARIOF R P，BEZERRA M O，VINZONOO S B．Dynamics of the saline front in the northern channel of the Amazon River:influence of fluvial flow and tidal range(Brazil)［J］．Journal of Coastal Research，2009，56:1414-1418．

［6］DOLGOPOLOVA E N．The conditions for tidal bore formation and its effect on the transport of saline water at river mouths［J］．Water Resources，2012，40(1):16-30．

［7］MENG Xia，PAUL M，CHRISTOPHER M，et al．Numerical simulation of salinity and dissolved oxygen at Perdido Bay and AdjacentCoastalOcean［J］．JournalCoastal Research，2011，27(1):73-86．

［8］HAN Zengcui，PAN Cunhong，YU Jiong，et al．Effect of large-scale reservoir and river regulation/reclamation on saltwater intrusion in Qiantang Estuary［J］．Science in China:Series B，2001，44(Sup1):221-229．

［9］史英标，潘存鸿，程文龙，等．钱塘江河口段盐水入侵的时空变化及预测模型［J］．水科学进展，2012，23(3):409-418．(SHIYingbiao，PAN Cunhong，CHENG Wenlong，et al．Temporal-spatial variation of salt water intrusion in the estuarine reach of Qiantang River and its numerical simulation［J］．Advances in Water Science，2012，23(3):409-418．(in Chinese))

［10］尤爱菊，韩曾萃，史英标，等．浙北引水对钱塘江河口段主要水环境的影响分析［J］．河海大学学报:自然科学版，2006，34(2):152-156．(YOU Aiju，HAN Zengcui，SHI Yingbiao，et al．Impact of water diversion from Xinanjiang Reservoir to north Zhejiang Province on water environment at Qiantang Estuary［J］．Journal of Hohai University:Natural Sciences，2006，34(2):152-156．(in Chinese))

［11］韩曾萃，程杭平．钱塘江江水含盐度计算的研究［J］．水利学 报，1981(6):46-50．(HAN Zengcui，CHENG Hangping．Study on the salinity calculation of the Qiantang Estuary［J］．Journal of Hydraulics Engineering，1981(6):46-50．(in Chinese))

［12］韩曾萃，程杭平，史英标，等．钱塘江河口咸水入侵长历时预测和对策［J］．水利学报，2012，43(2):232-240．(HAN Zengcui，CHENG Hangping，SHI Yingbiao，et al．Long-term predictions and countermeasures of saltwater intrusion in the Qiantang Estuary［J］．Journalof Hydraulics Engineering，2012，43(2):232-240．(in Chinese))

［13］曾剑，孙志林，潘存鸿，等．钱塘江河口径流长周期特性及其对河床变形的影响［J］．浙江大学学报:工学版，2010，44(8):1584-1588．(ZENG Jian，SUN Zhilin，PAN Cunhong，et al．Long-periodic feature of runoff and its effect on riverbed in Qiantang Estuary［J］．Journal of Zhejiang University:Engineering Science，2010，44(8):1584-1588．(in Chinese))

［14］潘存鸿，曾剑，唐子文，等．钱塘江河口泥沙特性及河床冲淤研究［J］．水利水运工程学报，2013(1):1-7．(PAN Cunhong，ZENG Jian，TANG Ziwen，et al．Sediment characteristics and study of riverbed erosion/deposition in Qiantang Estuary［J］．Hydro-Science and Engineering，2013(1):1-7．(in Chinese))

［15］潘存鸿，张舒羽，史英标，等．涌潮对钱塘江河口盐水入侵影响研究［J］．水利学报，2014，45(11):1301-1309．(PAN Cunhong，ZHANG Shuyu，SHI Yinbiao，et al．Study on saltwater intrusion in Qiantang Estuary affected by tidalbore［J］． JournalofHydraulic Engineering，2014，45(11):1301-1309．(in Chinese))

［16］潘存鸿，鲁海燕，曾剑．钱塘江涌潮特性及其数值模拟［J］．水 利 水 运 工 程 学 报，2008(2):1-9．(PAN Cunhong，LU Haiyan，ZENG Jian．Characteristic and numerical simulation of tidal bore in Qiantang River［J］．Hydro-Science and Engineering，2008(2):1-9．(in Chinese))

［17］TORO E F．Shock capturing methods for free surface shallow flows［M］．Chichester:John Wiley ＆ Sons，2001．

［18］ANASTASIOU K，CHAN CT．Solution of the 2D shallow water equationsusing the finite volume method on unstructured triangular meshes［J］．International Journal for Numerical Methods in Fluids，1997，24:1225-1245．

［19］于普兵，潘存鸿，谢亚力．二维风暴潮实时预报模型及其在杭州湾的应用［J］．水动力学研究与进展:A辑，2011，26(6):747-756．(YU Pubing，PAN Cunhong，XIE Yali．2-dimensional real time forecasting model for storm tides and its application in Hangzhou Bay［J］．Chinese Journal of Hydrodynamics，2011，26(6):747-756．(in Chinese))

［20］张泽锋，宋立松，汪荣胜，等．钱塘江河口盐度预报预警系统研制及应用［J］．浙江水利科技，2013(5):54-56．(ZHANG Zefeng，SONG Lisong，WANG Rongsheng，et al．Salinity forecast warning system in Qiantamg River estuary［J］．Zhejiang Hydrotechnics，2013(5):54-56．(in Chinese))