华东师范大学博士学位论文长江河口盐水入侵研究——北支倒灌、深水航道工程和冬季季 风的影响姓名：吴辉申请学位级别：博士专业：自然地理学指导教师：朱建荣20060401论文摘要本文采用三维数值模式和通量分析手段，定量研究了长江河口盐水入侵，特别是北支盐水倒灌及其动力机制，以及深水航道工程和冬季季风对盐水入侵的影响主要成果如下：( 1 ) 建立了一个高分辨率、移动潮滩边界的长江河口三维水动力、盐度数值模式运用模式对实测的潮位、流速、流向和盐度进行验证，验证结果良好 2 ) 对北支盐水倒灌进行了通量机制分析，结果表明L a 伊a 1 1 9 e 余流输运和潮泵输运是北支盐水进入南支并往下游输送的主要动力机制，在多数地区前者起主导作用，局部地区两者作用相当计算了南支各叉道的断面盐分通量以考查倒灌盐水团的输送途径，结果表明北支盐水出南北支分叉口后大部分沿白茆沙北水道下泻，少部分绕过白茆沙后沿着自茆沙南水道往下输送两部分盐水在七丫口附近河段汇合，然后其中一小部分通过南门水道进入新桥水道，并在新桥水道涨潮优势流的作用下向上游输移，其余盐水继续沿南支主槽下泻在新浏河沙以下，倒灌盐水通过南、北港分叉区域的几个叉道分别进入南港和北港。
在大潮期间，进入南港的盐量较多，在小潮期间，进入北港的盐量较多 3 > 采用北支上端的断面水体净通量作为倒灌强度的代表，在单日和半月两个时间尺度上研究了断面净通量和径流量、潮差的定量关系在单日时间尺度上，引入水体倒灌判别指标：细 青龙港潮差一1 - 5 8 ×e x p ( ' 2 4 5 7 ×! 鬻慕)加大于零时北支上端的断面净通量指向南支，倒灌发生；小于零时断面净通量指向北支，倒灌消亡；等于零时无净的水体通量基于如J D 得到了单曰断面净通量：单日断面净通量= 5 3 3 ．1 ×e x p ( o ．3 6 1 ×加D ) 一5 3 3 ．1 ∞s o ≥O )和单日断面净通量= 2 2 6 ．o ×如s 豇s o 锄加b )D l J > d m 帆D 1 ．1 J > d 聊加c )D ，J > d 砌拥，D I ．1 J d 玎d )Dz J 砌z 新，而且矢一l ，，一茧』> d 矗u ，j 点满足条件a ) 与b ) 或者a ) 与c ) 或者a ) 与d ) 的u 点为湿点，否则都为干点其中砌代表了干湿转化的水位梯度，本文中取5 c m ，在保证模式稳定的情况下越小越好。
同时满足a ) 和b ) 表示流速点以及其相邻的两个水位点的水深都大于临界值，故为湿点同时满足a ) 和c ) 或者同时满足a ) 和d ) 表示流速点及其一侧的水位点上的水深大于临界值，此时流速点的干湿由两侧水位梯度方向来决定当梯度指向干的水位点时流速点为湿点，否则就是干点对于v 点也一样，此处不再赘述对于一个给定的水位点，在D > d 埘胁时F s M 取1 ，否则为O 由上述判别法则我们可以看到，当流速点满足条件a ) 和d ) 时，水位点( i 0 ) 为干点，流速点( i j )为湿点，前者( F S M = 0 ) 的水深将由于后者导致的通量进入而发生变化，因此必需包括到计算区域内所以在数值模式中湿的水位点固然要参与到计算过程中，某些干点也必须添加到计算区域内一个简单的办法就是把所有的网格点不管干湿都计算一遍另一个方法同时也是本文所采用的就是引入一个新的开关变量F s M A D D 对于一个给定的水位点，如果它始终不可能参与到计算区域的时候取0 ( 如陆地) ，否则取1 这可以首先排除陆地上的网格点，减少了计算量2 ．1 ．6 数值求解方法模式中的动量方程分三步做时间积分：第一步： 些骘芋盟一璺≯一警一鱼磐吼叭Ⅳ砉( ≯盯南( 抄咄”∥出搿却衍2 一硝、^ ．’1 却、％”。
砚“ Ⅳ善c 彘，+ 量等舅p 筹拈一号} 毒j P ““a + n 以“( 2 - 4 3 )第二章长江河口三维盐水入侵数值模式的建立型堕；i 皇笪一堑坠笔笠一望堡笔竖一望! 乒一助砉( 幸一Ⅳ茜( 鲁卅潞矿万出髫却衍1 一并、^ ，’1 却、^ ，“’玎Ⅳ昙( 等) 一生巴罢7 宴；‰一壁也．晏- 一打+ D ‘峨n印h l h lp 铂；∞以蜘：’( 2 ·4 4 )第二步： 堕：云鳖= 寺砉k 等]△fD 街【d FJ 咝字丝= 古砉k 警]△fD ∞r ld J 』第三步：华呐驴笔△f一日f塑掣：也g D 掣△f—a ”( 2 —4 5 )( 2 - 4 6 )( 2 - 4 7 )( 2 - 4 8 )在第一步中，仅考虑非线性项、科氏力项、密度梯度力项、水平涡动糙滞项和由坐标变换产生的曲率项，所有这些项均作显式处理在第二步中，仅考虑垂向涡动粘滞项，并做隐式处理，采用三对角方法求解在第三步中，仅考虑外重力波产生的水位梯度力，这个快过程作隐式处理，把它代入连续方程，用半隐格式求解水位分布，消除了c F L 判据的严格限制，提高了计算效率求出水位分布后，由第三步可直接得出三维速度场的分布。
在数值求解盐度方程时，平流项的差分求解格式非常重要，常见的中央差格式会将下游的信息带往上游，从而造成物理失真和数值振荡，采用一阶迎风格式求解又会带来较大的数值耗散针对这个问题，本文采用欧拉一拉格朗日方法将这一方法引入到E C O M —s i 中的工作由朱首贤博士完成( z h ue ta 1 ．，2 0 0 1 ) ，下面简述这一方法的操作过程考察笛卡儿坐标系下的盐度输运方程：鲁+ “塞+ V 言+ w 龛= 差( 毛象) + c( 2 —4 9 )百+ “瓦+ ”瓦+ w 瓦2 瓦( ％夏J + ，s怛一4 w将f 2 —4 9 ) 左侧的平流项写成全微分的格式：粤= 要( 乜芸) + 只( 2 —5 0 )考察一流体微团在倒r 时刻的盐度s 矿，由( 2 —5 0 ) ，在∞+ 』Ⅻf 它移动到了新的位置，在移动过程中由于水平和垂向的扩散而使它的盐度发生了改变因此华东师范大学博士学位论文如果知道流体微团在上一个时刻的盐度，再计算它在运动过程中的垂向和水平扩散项就可以得到流体微团在当前时刻的盐度，即：品”1 = 品”+ ( 要( 毛要) + E ) ×&( 2 —5 1 )考察@ + 1 列} 时刻位于网格点( U 的的流体微团，则由( 2 —5 1 ) ，如果知道它在列f时刻的盐度，并将其垂向扩散项采用隐式计算，水平扩散项采用显式计算，就可以知道它在此时的盐度了。
流体微团在鲥f 时刻的盐度可以通过先确定它在剃f时刻的位置，再将蒯f 网格点上的盐度( 已求出) 插值到其上的方法求出整个过程可以用s 1 0 = P @ ”) + ( 要( K 要) ““+ E ”) x △r( 2 —5 2 )留芘来表示，s 蒜为铆+ ，Ⅻf 时刻网格点( 弛功上的盐度，s ”为列r 时刻的盐度场采用一次预估修正逆推得到网格点( “妨上流体微团在上一时刻的位置为：呻一厂_ ÷斗n + 1- ÷H + l1 ，’= ‘∽￡一o ．5 } ，( 1 小女一矿叫，女△f ) + 矿u ，女l △f( 2 ·5 3 ) Lj4n “呻 其中r t 为点( “的的空间位置，矿．* 为当前( “妨点的速度，一f 为时间步长，，为逆推得到的位置 ·) 为插值算子，把当前的流速或盐度值插到预估得到的位置上重复上一过程就是多次预估修正，但是采用一次预估修正已能获得较高的精度了将删f 时刻的盐度值插值到r ’上得：s 观= 妒( s ”) + ( 罢( 蚝象广1 + c ”) ×心= ，( ≯) + ( 丢( 邑象厂1 + E ”) ×出( 2 - 5 4 )插值算子，( ·) 对计算的精度十分重要，本文在水平方向上采用一阶、二阶混合的拉格朗日插值方法，垂向上采用线性插值方法。
这种方法即提高了精度，又有效的避免了二阶拉格朗日插值引起的附加数值振荡，具体参考( z h us h o u X i 趴e t a l ．，2 0 0 1 ) 第二章长江河口三维盐水入侵数值模式的建立2 ．2 模式的设置2 ．2 ．1 模式的范围和网格图2 ．1 模式的计算网格，图中大范围计算网格分辨率为实际的三分之一．2 3 ．华东师范大学博士学位论文模式计算范围包括整个长江口、杭州湾，东边到1 2 4 ．5 0 E 附近，北边到3 3 0 N附近，南边到2 8 0 N 附近在模式中各地理要素的坐标一律转化为1 9 5 4 北京坐标河口地区水流为双向流，而当上游边界以径流量给出的时候只能提供单向流，因此上游边界应该尽可能的取得远一些理想的状态是将上边界取在大通站，这样就不存在单向流和潮汐振荡的影响，但这样一是增加计算量，二是对研究区域的意义不大，本文中将上游边界延伸到江阴以上1 5 0 k m 南北支分叉口区域的网格加密，在北支上端的分辨率为7 5 m 左右；深水航道区域网格加密且拟合导堤、丁坝，航道区域的分辨率为1 0 0 m 左右；南槽沿岸附近网格加密，网格分辨率为7 0 m 左右，口内其余部分分辨率在1 0 0 m 至5 0 0 m 不等，口外较疏，最大为7 k m左右。
垂向采用仃坐标，均匀的分成5 层2 ．2 ．2 长江口及临近海区的地形2 ，u) ⋯u⋯u⋯) ub I ^ J n l x T a ⋯r k I n 、图2 ．2 长江口、杭州湾及邻近海区水深本文在长江口内主要采用2 0 0 3 年长江河口水文观测得到的水深资料，其中局部地区采用更新的资料长江口外的水深变化不甚剧烈，本文沿用以前海图数字化的资料( 图2 ．2 ) 所有水深资料都校正到8 5 基面姗舢删姗姗姗姗珊瑚善!^§，i3第二章长江河口三维盐水入侵数值模式的建立2 ．3 ．3 外海开边界条件外海开边界水位采用8 个天文分潮M 2 ，S 2 ，N 2 ，K 2 ，K 1 ，O l ，P l 和Q l 驱动8 f = 氢+ ∑工日，c o s ( 国，f + ( 矿．+ “) 一g ，)( 2 —5 5 )』2 l其中f 为潮位，氏为余水位，厂为节点因子，珏为振幅，国为角频率，昏为迟角K + 挑为订正角，它可由具体的年、月、日求得图2 ．3 大模式计算网格为了提供合理的开边界条件，建立了一个包括整个渤、黄、东海的大模型来模拟各个天文分潮，模式范围西至1 1 7 ，5o E ，东至1 3 1 0 E ，南至2 3 ．5 0 N ，北至4 1 0 N( 图2 ．3 ) 。
大模式较成功地模拟出了8 个分潮在东海黄海渤海的转播特征( 图2 ．4 、2 ．5 ) ，再现了计算区域内半日分潮的5 个无潮点和2 个蜕化的半个无潮点、全同分潮的3 个无潮点与6 5 个潮位站的观测资料比较，模式计算的M 2 、s 2 、华东师范大学博士学位论文K 1 和0 1 分潮振幅和位相的均方差分别为7 ．8 5 c m 和6 ，8 1 ( 表2 ．1 ) 、5 ．0 4 c m 和8 ．1 4 、3 ，8 2 c m 和1 3 ．0 4 、4 ．3 4 c m 和9 ．3 3 与1 7 个潮位站的观测资料比较，N 2 、K 2 、P 1 和Q 1 分潮振幅和位相的均方差分别为3 6 4 c m 和7 ．8 9 、8 ．4 7 c m 和1 0 ．5 1 、1 ．7 6 c m 和7 ．5 6 、1 ．5 0 c m 和2 6 ．3 4 模式模拟的结果可为本文小区域模式提供可靠的外海开边界潮汐调和常数资料⋯⋯“2 ”3 焉”⋯⋯⋯⋯⋯””⋯⋯⋯2 ⋯2 ⋯3 啬．∞⋯⋯⋯⋯”1”⋯⋯⋯⋯”j 高⋯”⋯⋯⋯11 ””⋯⋯⋯2 ⋯j 篙”6 ”7 ”⋯⋯⋯图2 ．4 大模式计算的半日潮同潮图，实线振幅，虚线位相，左上、右上、左下和右下分别代表M 2 、S 2 、K 2 和N 2第二章长江河口三维盐水入侵数值模式的建立⋯⋯⋯⋯”乙啬。
箍⋯⋯“⋯o⋯⋯⋯⋯2 ⋯。