钱塘江河口洪水特性及动床数值预报模型
摘要:在分析钱塘江河口水文、地形及洪水位实测资料的基础上,得到了钱塘江河口段潮汐、径流、河床三者的非线性关系对河口段的洪水位有显著影响的,以此为基础建立了洪水预报的一维动床数值预报模型。模型复演了钱塘江河口一次典型的洪水位过程,并利用动床模型和定床模型的两种模式对各站洪水位进行了比较,两种模式计算得到的沿程洪水位可差0.3~2.0m,但动床模拟与实测过程十分吻合,表明在象钱塘江河口这样大冲大淤的洪水位预报中考虑动床冲淤过程是必要的。
关键词:钱塘江河口 洪水特性 动床模型
0前 言
钱塘江河口是一个冲积性河口(图1),洪冲潮淤,年内基本平衡。每年汛期前4月份,如前期为枯水,江道淤高,汛期即使遭遇中小洪水,沿程洪水位也比常年高出甚多,即小流量高水位的现象在钱塘江河口非常多见,这一特点不同于其它河流,它给洪水预报提出了极高的要求。对于河床易冲易淤的钱塘江河口,洪水位预报的难点主要有二,第一是由于钱塘江河口河床地形受径流丰枯多变,洪水预报前期地形数据难以准确获得,一般每年4月、7月、11月份施测三次地形,而洪水一般多发生在6月中旬至7月中旬的一个月内,如4月至6月中旬上游雨量偏多或偏少,钱塘江河口就会发生较大幅度的冲淤变化,如仍采用4月份的江道地形预报6月中旬以后发生的洪水过程,会带来很大的误差;第二是必须考虑洪水过程中河床的冲淤变化,原因在于钱塘江河口由细粉沙组成,床沙中值粒径在0.02~0.04mm,泥沙起动流速0.5~0.7m/s,极易起动,致使河床在洪水期变化较快,一般历时3~5天的洪水可将前期7个月淤积的泥沙搬运至河口段的下游。这表明在钱塘江河口洪水预报中必须考虑河床的冲淤过程。
图1钱塘江河口示意图
Fig.1 Sketch map of Qiantang Estuary
洪水预报的主要任务是根据河段上游的已知流量过程,预测下游沿程重要城镇的流量、流速及水位过程。洪水系长波,其运动规律可用圣维南方程组来描述。钱塘江河口段河床由细粉沙组成,易冲易淤,河床惯性较小,特别是在洪水期较短的时间内,河床冲刷剧烈,同一洪峰流量的洪水位可相差1m之多,这主要是由于在钱塘江河口段潮汐、径流、河床地形三者的非线性关系对河口段的洪水位有显著的影响。几十年来,工作者在分析钱塘江河口大量水文、潮汐及地形的实测资料的基础上,得到了影响河口洪水位的几种因素,并根据河口段潮流、泥沙运动与河床变形理论进行了数值预报模型的研究,取得了不少的研究成果[1][2]。但限于当时计算技术条件,对洪水前期的江道地形采用经验修正的方法得到,再用经验合轴相关曲线作闸口站洪水位预报,它的优点是方便快速,缺点是不能动态预报洪水位沿时空的变化过程。本文在前人工作的基础上,建立了一维动床数值预报模型,进行了钱塘江河口的一次典型洪水过程的数值预报,其速度和精度可满足预报要求,为实时预报奠定了基础。
1钱塘江河口洪水位变化特性
1.1钱塘江河口洪水特性
钱塘江河口七堡以上河段最高水位受梅汛期洪水控制(仅个别年份由台风暴潮引起)。由于钱塘江河口具有潮大流急,河床冲淤幅度大、速度快的特点,使得该河段的洪水位变化规律既不同于无潮河流,也不同于冲淤变化较小的潮汐河流。具体地说,它除了与汛期洪水流量大小直接有关外,还受江道冲淤面貌、下游潮汐、河口下段尖山河湾主槽走向、长度及人类活动(包括新安江水库的兴建及治江围垦)等诸多因素控制。一般地,闸口洪水位(Z)可由闸口洪峰流量(Q)、闸口至盐官江道容积(V)、下游澉浦高潮位(ZZ)三者的非线性回归关系确定,由此建立经验预报合轴相关图[2]。
1.2 泥沙输移对洪水位影响
钱塘江河口泥沙输移对沿程洪水位影响比较敏感。实测资料分析表明,在枯水大潮期间,盐官以下河段落潮含沙量大于涨潮含沙量,而盐官以上河段涨潮含沙量大于落潮,净泥沙输移是指向上游,从而造成钱塘江河口自上而下的淤积,导致该河段年内潮淤洪冲,年际间表现为丰水年冲刷,枯水年淤积的格局。泥沙输移的最终结果使闸口至盐官河段河床发生冲淤变化,进而影响行洪面积。该河段的冲淤面貌可用钱塘江河口沙坎高程及闸口至盐官河段7m下容积等指标表征,它们均可用闸口站低水位集中反映。从表1的几次洪水资料可以看出江道冲淤面貌对洪水位的影响起决定作用。比较1955年6月与1968年7月两次闸口洪水位,由于汛前江道面貌不同,1955年汛前江道冲刷,面积、容积均很大,汛前闸口断面低水位仅3.87m,沙坎高程2.1m,而1968年汛前遭遇枯水而淤积,江道面积、容积较小,仅2.68亿m3,闸口低水位高达6.14m,沙坎高程4.0m,尽管闸口洪水流量相差2倍多,但其最高水位非常接近,分别为9.07m和9.01m,出现小流量高水位的现象主要在于江道淤积所致。而1968年7月与1955年4月另一次洪水相比闸口洪峰流量相近,分别为12000m3/s和11000m3/s,但由于江道面貌不同,两者的洪水位差达2m之多。
表1 实测典型洪水比较
Tab.1 Comparison of measured typical floods in the Qiantang river
年.月.日 | 洪峰 流量/m3.s-1 | 闸口 洪水位/m | 汛前闸口低水位/m | 河口沙坎 高程/m | 澉浦 高潮位/m | 闸口至盐官 容积/亿m3 |
1955.6.22 | 31000 | 9.07 | 3.87 | 2.10 | 6.61 | 9.31 |
1968.7.10 | 12000 | 9.01 | 6.14 | 4.0 | 6.60 | 2.68 |
1955.4.17 | 11000 | 7.00 | 4.47 | 3.70 | 4.56 | 7.98 |
1997.7.9 | 15000 | 9.67 | 5.88 | 4.5 | 6.31 | 3.05 |
1.3 潮汐对河口洪水位的影响
在钱塘江河口洪水演进中,下游澉浦潮汐大小对该河段洪水位起顶托作用,相同的江道条件与洪水条件下,潮汐大则洪水位高。如1994年6月18日与1995年6月25日两次洪水,上游富春江电站的洪峰流量基本相同,分别为10100m3/s和10600 m3/s,洪水前闸口至盐官7m下容积也基本相同,均为4.0亿m3,但下游澉浦的高潮位不同,分别为5.15m、4.31m,差0.84m,而闸口的洪水位分别为8.1m、7.56m,两者相差0.54m之多。可见潮汐对洪水位影响比较明显。
2 动床数值预报模型的建立
2.1 基本控制方程及定解条件
钱塘江河口经治江围垦后,河道已相对比较规则,可以用一维动床数学模型来描述水流运动,考虑到水体挟带的泥沙主要以细颗粒悬移质为主,河床冲淤也以悬移质起主要作用,因此本文在模型中也主要考虑悬移质泥沙,单一河道悬沙数学模型基本控制方程可以表述如下:
水流连续方程 (1)
水流运动方程: (2)
悬沙输移方程: (3)
河床变形方程: (4)
挟沙能力公式: (5)
式中:x为河道里程,t为时间,Q为流量,B为河宽,Z为潮位,A为断面过水面积,为谢才系数,R为断面水力半径,S为断面平均含沙量,S*为断面平均挟沙能力,、为泥沙冲淤系数,ω为泥沙沉降速度,为断面平均河底高程,为单位河段傍侧入流量。E为泥沙纵向扩散系数,γs为泥沙干容重。
将方程(1)~(3)式按Preissmann四点偏心隐格式离散,方程(4)用一般的时间差分格式离散得到差分方程。限于篇幅,基本方程的定解条件及离散方程详见[3][4]。
但象钱塘江之类的河口,在径流和潮流作用下形成往复水流,就流动方向而言,存在如图2的四种情况。(a)、(b)两种情形与单向流动没有什么区别,水流和含沙量方法也与单向河道相同;(c)、(d)两种情形与单向流动的差异在于水流和含沙量内边界提法不同,特别是(d)情形,计算河段的上、下边界不需要含沙量边界条件,其计算边界条件应在停滞点(Q=0)所在断面给定水流和含沙量的内边界条件。对于(c)不需给含沙量内边界条件,只要将河段分成两单向河段处理即可。下面重点说明(d)情况的水流、含沙量内边界条件如何给定。
图2 潮流往复流动的四种情况
Fig.2 Sketch of the tidal flow
在水流运动方程中的阻力项已作了单向流动的假定,正、负两种方向的流动阻力系数不等,在包含正、负二向流动的河段中,应将二向流动河段分别处理。假定滞流点(Qk=0)在Δt内位置固定不变,m~k和k~n两段分别可用Preissmann格式离散得离散方程,其式与在m~n断面直接离散得到的离散方程的形式区别仅在于用△Xmn/2代替其中的。根据Qn,Qm流量的大小,利用线性插值确定滞流点的位置:
(6)
滞流点的含沙量可由下式确定:
(7)
经离散得到:
(8)
这样以K断面为界,可按两股单向流动分别计算沿程断面的含沙量。因此计算内边界点含沙量时必须首先求得内边界点的位置和数量,然后(8)分别确定内边界点的含沙量。
2.2 关键问题的处理
动床数学模型计算的关键是准确选用水流、泥沙基本方程的有关参数。
①床面阻力系数
河床糙率沿程分布采用实测水位资料率定得到,涨、落潮糙率取不同值。闻家堰以上河段涨、落潮糙率取值范围为0.02~0.040;闸口以下河段涨潮糙率为0.01,落潮糙率0.01~0.012;闻家堰至闸口河段涨、落潮糙率取值范围为0.01~0.02。
②水流挟沙能力
目前,挟沙能力仍以据本河段实测资料得到的经验关系最为可靠,通过大小潮及洪水实测资料分析得到钱塘江河口段水流挟沙能力公式为[1]:
(9)
式中m=0.95~1.1,H为水深,K为挟沙能力系数。K值的合理选定可从实测水沙资料推求,考虑到沿程水流、泥沙等条件的不同,其值也有所不同,而且在丰枯及潮水期也有一定的差别。因此在具体选择K值时可通过含沙量及河床变形率定和验证得到。
③泥沙冲淤系数T1、T2
本文中的泥沙冲淤系数系底部含沙量与断面平均含沙量、底部挟沙力与断面平均挟沙力的比值,根据林秉南、韩曾萃的研究成果[5]:
, (10)
④河床冲淤判别条件
采用含沙量与挟沙力对比的判别条件,即当S>S*时河床淤积,当S<S*时,且U>Uc,河床冲刷。式中:Uc为泥沙起动流速,可用武汉电力大学的起动流速公式计算[6]。
2.3冲淤面积的分配
一维输沙模型给出了断面的冲淤面积和河段的冲淤量,要进一步了解冲淤沿纵向和横向的变化,冲淤量的合理分配是保证一维动床数学模型精度的根本问题。冲淤量的横向分配最常用的是根据水沙条件和河床组成的某些数学关系对河床进行分配,主要有平均分配法、按面积比考虑分配、按能量比进行分配、按切应力大小进行分配、按水流挟沙的饱和程度考虑等[7]。本文采用冲淤面积沿水深分配的方法,即:△Zsi=hij△Asi/Ai,,这种分配方法的主要缺点是不区分滩槽且不考虑水沙条件及床面切应力条件。在潮汐河口河床冲淤计算中冲淤量的纵向分配也十分重要,由于在长河段长时段的演变计算中上、下河段水动力和泥沙条件的不同,其冲淤性质可能出现如下几种情况:
上、下游河段冲淤性质相反——上段冲刷、下段淤积和上段淤积、下段冲刷;
上、下游河段冲淤性质相同——上下河段同为冲刷和上下段同为淤积;
对上述几种情况要分别处理。如上、下游河段冲淤性质相反,则首先要确定冲淤量为0的断面位置,求得该点的含沙量作为上游断面的出口含沙量和下游河段的进口含沙量,分别计算冲淤量,如上下河段冲淤性质相同,情况较为简单,可从河床演变自动调整机理出发,断面冲淤量取决与断面含沙量与挟沙力的差,利用上述沿程冲淤量的处理方法即可反映钱塘江河口沙坎的纵向演变过程。
2.4 区间入流概化
钱塘江河口段的洪水不仅包括富春江电站的下泄洪水流量,而且还包括沿江各支流(主要有分水江、浦阳江、壶源江、禄渚江等)汇入干流的洪水流量。分水江、浦阳江、壶源江、禄渚江的控制站分别为五里亭、诸暨、高峰及徐畈。因此本预报模型将壶源江、禄渚江的洪水及其它沿江分散集雨面积的洪水按其控制面积的大小分摊到分水江和浦阳江。在此基础上,利用动床洪水预报模型联合求解五里亭~桐庐河段、富春江电站~桐庐以及桐庐~澉浦河段的洪水演进过程,其中桐庐节点为分水江和富春江的汇合点,浦阳江口节点为浦阳江和富春江的汇合点,这些节点满足水量守恒、能量守恒和沙量守恒条件。
3 模型验证
一个完整的动床冲淤数学模型应包括枯水期大、中、小潮淤积过程以及洪水的冲刷过程的验证计算。其中淤积过程的验证计算见[8],下面重点介绍洪水短期的冲刷过程。1997年7月7日至12日,钱塘江河口发生了一次典型的洪水,由于汛前1996年11月至1997年4月,富春江电站下泄流量平均为442 m3/s,比常年减少30%,致使河口江道淤积,容积较小,汛前4月份江道容积仅3.13亿m3,又5月~6月径流仍较小,江道持续淤积,估计在“7.9”洪水前江道容积约2.8亿m3,致使1997年“7.9”洪水闸口流量仅15000 m3/s,约五年一遇,钱塘江河口沿程洪水位超记录,且高水位持高不下的危急局面。本模型以此为例复演了钱塘江河口沿程洪水位过程。
3.1 洪水位验证结果
动床模型的计算结果见图3。由图3可知:由动床计算得到的各站水位过程与实测过程十分吻合,闻堰以下河段最大误差在15cm内。在此基础上再用定床模型进行了计算,两者的结果比较见表2。
表2 定床和动床洪水计算结果比较(m)
Tab.2 Comparison of calculated results between fix-bed and mobile-bed model
水位站名 | 桐 庐 | 富 阳 | 闻 堰 | 闸 口 | 七 堡 | 仓 前 | 盐 官 | |
高水位 | 动床 | 16.17 | 11.74 | 10.20 | 9.82 | 9.12 | 8.57 | 7.64 |
定床 | 16.48 | 12.96 | 11.82 | 11.52 | 10.13 | 8.74 | 7.77 | |
差值 | 0.31 | 1.22 | 1.62 | 1.70 | 2.01 | 0.17 | 0.13 | |
洪水末 | 动床 | 12.85 | 9.12 | 8.11 | 7.72 | 7.11 | 6.28 | 4.33 |
定床 | 13.46 | 11.08 | 10.5 | 10.31 | 9.17 | 7.35 | 4.81 | |
差值 | 0.61 | 1.96 | 2.39 | 2.61 | 2.06 | 1.07 | 0.48 |
从图、表可知动床与定床计算得到的水位过程各站差别很大。其中,桐庐站最高水位动床和定床计算差别不大,仅为0.31m,富阳~七堡河段高水位可差1.22~2.01m,七堡以下河段仓前~盐官河段高水位主要受下游高潮位控制,动、定床模型的差异较小,约0.13~0.17m,但洪水后期动床和定床差别较为明显。其原因在于该河段的洪水位主要取决于潮汐、径流、地形的非线性作用,特别是洪水期河床的冲刷作用对其水位的变化过程有显著的影响,其中上游桐庐段,由于河床为粗沙组成,动、定床计算对高水位的影响较小,说明可以不必考虑动床;闻堰、闸口、七堡及仓前动床和定床的计算对洪水位过程相差很大,主要是钱塘江河口河床惯性小,易冲易淤,显然动床模型更符合实际,说明本河段必须用动床计算;而富阳站受下游动床计算的影响而变动;仓前以下的高水位则受下边界澉浦高潮位控制,动床与定床相差不大,而低水位则又受河床冲淤影响差别较大。
3.2 泥沙冲淤时空分布
表3是河段冲淤量实测与计算的对比,由表可知实测与计算基本吻合。图4是河口段沿程各断面河床高程变形速度(dz/dt)与河段冲淤量随时间的变化过程。
表3 1997年“7.9”洪水冲刷量验证(亿m3)
Tab.3 Calibration of amount of erosion
河 段 名 称 | 1997年4月 | 1997年7月容积 | 冲 刷 量 | ||
实 测 | 计算 | 实测 | 计算 | ||
闸口至七堡 | 0.68 | 1.37 | 1.26 | 0.69 | 0.58 |
七堡至仓前 | 0.65 | 1.17 | 0.99 | 0.52 | 0.34 |
仓前至盐官 | 1.64 | 1.83 | 1.87 | 0.19 | 0.23 |
从图4可知在洪水作用下闸口至七堡河段河床为单向冲刷,冲刷强度随洪水流量的减小及断面的冲深而减弱,洪水流量最大时,河床冲刷速度亦最大,河段单位时间的冲刷量洪水涨水期大于落水期;盐官断面则在洪水与潮水的双重作用下河床时冲时淤,表现为潮淤洪冲,但由于洪水占主要地位,最终结果为河床冲刷;盐官以下断面由于潮水占主要地位,最终结果为河床淤积。
图3 沿程各站水位实测与计算值对比图
Fig.3 Comparison between the Calculated and the measured water level
图4 冲淤沿时空分布
Fig.4 The space and time distribution for amount of deposition/erosion
4 动床数值预报的讨论
影响洪水成果的主要因素有:洪水流量(包括支流入汇)、河床初始地形及模型中的有关物理参数,前面二个涉及模型的初边值条件的确定,可用流域水文学的降雨径流模型求得。江道初始地形可选洪水较近的实测地形资料为基础尔后用一维动床数学模型修正得到。下面重点讨论模型的物理参数,即水流挟沙能力(系数K)、泥沙冲淤系数(T1,T2)及泥沙沉降速度(ω)对洪水位计算成果的敏感性。
潮流挟沙能力:模型采用的形式采用式(9),对其系数K=6,12,18三种情况进行了对比计算,结果表明系数K对闸口的洪水位影响较大。其中K=12、18计算过程与实测过程比较吻合。K值太小,洪水的冲刷能力小,洪水位高,持续时间长。经钱塘江河口多种地形下的多次洪水位预测成果并与实测值的比较得到k值的一般取值范围为当闸口至盐官河床容积大于5亿m3,可采用定床模型计算,定床与动床差异较小;容积在3.8~5亿m3,挟沙力系数采用k=6~8;容积在3.8亿m3以下,系数k=12~14。主要原因在于水流挟沙力系数与床面泥沙沉积特性有关。
泥沙冲淤系数(T1,T2):泥沙冲淤系数(T1,T2)是反映床沙与悬沙交换使挟沙水流从非饱和向饱和恢复速度的一个综合系数,本文对冲淤系数(T1,T2)=0.2,1.0及>1.0进行了比较。当(T1,T2)<1.0时,闸口计算洪水位与实测过程相差很大;(T1,T2)>1时,计算洪水位和含沙量与实测过程非常吻合。
泥沙沉降速度(ω):钱塘江河口悬沙沉降速度受江水含盐度、含沙量等影响外,在动床冲刷计算中还受床沙粗化等因素影响。因此,本文比较了泥沙沉降速度为ω及1.5ω时闸口的洪水位,其中ω为初始时刻的各断面泥沙沉降速度。计算结果表明泥沙沉降速度对洪水位影响较小,两者的差异平均值仅在7cm左右。因此,与水流挟沙能力系数及泥沙冲淤系数相比,泥沙沉降速度的变化对钱塘江河口洪水位预报误差的影响较小,其沉降速度可按床沙中值粒径的实际值给定,而不必考虑冲淤过程中泥沙粒径变化引起的沉降速度的变化。
5 结 论
通过实测资料得到钱塘江河口段潮汐、径流、河床地形三者的非线性关系对河口段的洪水位产生显著影响的,同一洪峰流量和潮汐条件下河口段洪水位可差1~2m,表明在河床易冲易淤的钱塘江河口洪水位预报模型中考虑河床冲淤演变的必要性。在此基础上建立了洪水位动床数值预报模型,模型复演了钱塘江河口1997年“7.9”洪水过程,计算得到的洪水过程、冲刷量及时空分布与实测值吻合较好。最后,通过数值试验对预报模型的敏感性进行了分析,水流挟沙能力(系数K)对洪水位计算的敏感性较大,泥沙冲淤系数(T1,T2)的大小其次,而泥沙沉降速度(ω)对河口段沿程洪水位的预报影响较小,在钱塘江河口洪水预报中可不必考虑河床冲淤过程中泥沙粒径变化引起的沉降速度的变化。
[1]韩曾萃,程杭平.钱塘江河口河床变形计算方法及其应用[J].泥沙研究,1987,(3): 43-54.
[2]韩曾萃,周文波. 钱塘江河口段最高水位的研究[R]. 浙江省河口海岸研究所技术报告,1987.2.
[3]Cunge.J.A.Holley.F.M.JR. and Verywey.A.(1980) Practical Aspects of Computation River Hydraulics[M]. Pitman,London U.K.
[4]史英标,韩曾萃. 动床模型在钱塘江河口洪水预报中的应用[C]. 河流模拟的理论与实践,武汉电力大学出版社,1998.10.
[5]林秉南,黄菊卿,等. 钱塘江河口潮流输沙数学模型[J]. 泥沙研究,1981,(2):16-17.
[6]谢鉴衡,等. 河流模拟[M],电力出版社,1993.10.
[7]杨国录. 河流数学模型[M]. 海洋出版社,1992.
[8]史英标. 钱塘江河口一维动床模型研究[R]. 浙江省河口研究院内部报告,1997.1.