长江口径潮相互作用下洪枯季余水位变化

翟 玮,张 蔚,季小梅,陈 婷

(1.河海大学 江苏省海岸海洋资源开发与环境安全重点实验室,南京 210024; 2.上海滩涂海岸工程技术研究中心,上海 200061; 3.河海大学 疏浚技术教育部工程研究中心,南京 210024)

1 研究背景

河口是陆海交汇的过渡地带,受到流域径流和海洋潮汐的共同作用,水动力条件复杂。余水位即潮周期内的平均水位,滤去了周期性变化的潮汐,其时空变化反映了径潮动力间的非线性作用[1]。根据一维动量守恒方程,余水位梯度项基本与摩阻项相平衡[2]。为探究径流、潮汐及其相互作用对余水位的影响,多从动量方程的摩阻项入手,将其分解为径流作用、径潮相互作用和潮不对称作用产生的摩阻项[3-6]。近年来随着技术手段的发展,学者们提出了新的方法来研究径潮间的非线性作用,包括连续小波变换[7-8]、短期调和分析[9]、非稳态调和分析[10]等。这些分析方法结合摩阻项的分解,有助于解释评价径流、潮汐及其相互作用对余水位的相对重要性。如Buschman等[4]基于连续小波变换方法和对摩阻项中非线性流速乘积的一阶及三阶展开,建立线性回归模型,通过潮平均摩阻的变化预测余水位的变化。

长江口是世界第三大河口,径流量巨大且季节性变化特征明显,潮汐强度中等,径潮相互作用显著,深刻影响着此处的余水位。过去对长江口余水位的研究,大多采用实测资料的统计分析方法。刘晓婉等[11]发现崇明岛南部近岸水位主要由潮汐驱动形成。郭磊城等[12]基于实测径流和潮位资料发现,长江口洪季径流量的增大使得余水位显著抬高,洪枯季差异从上游大通到下游牛皮礁逐渐减小。杨正东等[13]发现长江口余水位的变化主要与径流量有关。此外,宋永港等[14]基于数值模拟,对长江口北支潮位潮差时空变化展开研究,发现1—7月份的月平均潮位逐渐增大,而8—12月份逐渐减小,该变化趋势由径流产生的余水位决定。综上而言,针对径潮动力下余水位时空变化的研究仍有待深入。余水位影响着水位与水深,同时作为径潮动力非线性的典型结果,在长江河口动力格局变化的当下也成为分析径潮相互作用的有效切入点,对防洪安全、水位预报、水资源管理等科学问题具有重要意义。

本文基于二维水动力数值模型,获取长江口主要站点高频、长期的潮位和流速数据,分析余水位时空变化特征;利用连续小波变换分解流速信号,通过摩阻公式计算潮平均摩阻的各分量,结合简单线性回归模型分析径潮动力相互作用下余水位的变化机制。

2 研究区域及研究方法

2.1 研究区域

长江口是衔接了长江和东中国海的主要通道,全长约660 km,上至安徽大通,下至外海50 m等深线附近。长江口整体平面形态呈现为“三级分汊、四口入海”的复杂格局,河口呈典型喇叭型,上游宽度13 km,外海口门处宽度90 km,其地形及主要站点位置如图1所示。由于属亚热带季风气候,河口径流季节性变化显著,5—10月份径流量占年径流量71.7%,11月份至次年4月份径流占28.3%[15]。全日分潮和半日分潮是进入河口的主要分潮,其振幅比AD1/AD2约为0.24,河口内潮汐类型属于非正规半日潮。

图1 长江口地形及主要站点位置Fig.1 Map of the Yangtze River Estuary and locations of hydrologic stations

2.2 研究方法

2.2.1 水动力模型

本文采用课题组已建立的长江口整体平面二维水动力数学模型,对径潮动力相互作用下余水位时空变化展开研究。模型采用非正交的三角形网格,覆盖整个长江口区域,陆地边界能较好地贴合实际岸线,模型网格总数为157 311,网格节点数为75 154,网格大小在50～10 000 m之间变化,精度满足需要。模型地形资料采用长江水利委员会长江下游水文水资源勘测局测得的长江口水下实测地形,并通过14个潮位站点和5个流速站点的实测数据进行率定与验证。大多数测站水位和流速数据误差均在10%以内,结果良好,因此可认为模型模拟结果能较好地反映长江口的水动力条件。具体率定与验证过程见文献[16]。本文模型参数设置如表1所示。

表1 模型参数设置Table 1 Setting of model parameters

为了充分消除潮汐的影响,模型上边界设置在枯季的潮区界大通,并将2016年大通实测径流量作为上边界条件;外海边界则设置在径流作用完全消失的124.5°E处,分潮潮位由东中国海潮波数学模型[17]提供。模型采用冷启动方式,对2016年长江口的潮位和流速开展模拟。

2.2.2 连续小波变换及潮波分解

作为非稳态调和分析的重点方法之一,连续小波变换(Continuous Wavelet Transform,CWT)通过对母小波进行缩放和平移,对潮波信号进行了多尺度的细化分析,利用小波变换在一个带通滤波器上对潮波信号进行处理,得到信号随着时间和频率的变化过程[18-19]。

CWT方法的特点在于它可以比较准确地得到不同周期分潮簇的振幅与相位,非常适应潮汐频带非几何变化间隔的特性;其缺点是不能较好地区分同一分潮簇内的不同分潮[8],且对小波的选择具有较强的依赖性。

在径潮动力相互作用的背景下,CWT方法可以显露潮波对变化的径流量近乎瞬时的响应[9]。利用CWT方法处理流速数据,可将原始流速时间序列分解为不同时频域信息,将同一周期的分潮簇结果归类合并,从而得到主要的全日潮簇D1、半日潮簇D2以及四分之一日潮簇D4流速的振幅(U1、U2及U4)和初相位(φ1、φ2及φ4)。

2.2.3 摩阻效应

根据Buschman等[4]提出的针对余水位变化的新的分析方法,可将动量方程中的摩阻项按照其来源分解为径流作用、径潮相互作用和潮不对称作用产生的3个不同分量Sr、Srt及St。该方法的核心在于利用切比雪夫多项式对摩阻项中的非线性乘积U|U|进行一阶和三阶展开,并将流速时间序列分解后经无量纲化处理[5]代入,即:

(1)

(2)

对摩阻项在一个全日潮周期内进行积分,整理后得到公式如下:

(3)

(4)

(5)

日均流速和分潮簇流速振幅在经过无量纲化处理后均<1,结合式(3)—式(5)分析可知,Sr、Srt和St大致范围在0～1范围内,反映各项作用对潮平均摩阻的贡献。其中,Sr只与日均流速有关,因而表达的是径流作用对潮平均摩阻的贡献;Srt与平均流速和各分潮簇流速振幅的乘积相关,反映径潮相互作用对潮平均摩阻的贡献;而St则与分潮簇流速振幅和相位有关,反映潮不对称作用对潮平均摩阻的贡献。

3 余水位时空变化特征

根据二维水动力模型模拟结果,对长江口主要站点的模拟水位取日均值,从而滤去全日潮簇及更高频分潮簇潮位对水位周期性的影响,得到余水位。图2为2016年大通站的日均径流量及主要站点的余水位。

图2 2016年大通站日均径流量及主要站点余水位Fig.2 Average daily runoff at Datong station and residual water level at major hydrological stations in 2016

从图2(a)可以看出,长江口的径流量季节性变化明显,洪季(5—10月份)流量明显大于枯季(1—4月份及11—12月份)流量,日均径流量最大可达70 000 m3/s。根据图2(b)可知,长江口余水位从下游吴淞向上游南京逐渐抬升,具有稳定的梯度。在枯季,余水位沿程增幅较小,以两周为周期的大小潮变化始终明显,大潮时余水位要高于小潮时;在洪季,余水位沿程增幅较大,在上游站点其变化与径流量的变化十分相似,往下游余水位大小潮变化逐渐明显。上述长江口余水位的时空变化体现了径潮动力的共同作用,但其变化原因还需通过理论方法作出解释,具体见本文5.2节末尾。

4 分潮簇流速振幅时空变化

基于二维水动力模型的模拟结果和连续小波变换方法,对各站点流速时间序列进行调和分析,得到全日潮簇(D1)、半日潮簇(D2)和四分之一日潮簇(D4)流速的调和常数。图3为各分潮簇流速振幅的洪枯季沿程分布。

图3 各周期分潮簇流速振幅洪枯季沿程分布Fig.3 Distribution of velocity amplitude of tidal clusters in different periods during wet season and dry season

由图3可知,长江口内半日潮簇D2流速振幅最大,枯季在吴淞处可达0.9 m/s,全日潮簇D1和四分之一日潮簇D4流速振幅大致相当,均<0.2 m/s。潮波从下游的吴淞向上游传播过程中,受径流顶托和摩阻耗散,各分潮簇流速振幅大体上沿程衰减。其中在吴淞至徐六泾段,断面收缩明显,潮能集中,因此流速振幅衰减较慢。与其他两项不同,D4为浅水分潮,其能量也能通过浅水非线性作用从天文分潮转移而来,因此在吴淞至徐六泾段流速振幅明显增大。

对比洪枯季分潮簇流速振幅,各站点枯季流速振幅均大于洪季,可以得知洪季更大的径流量加剧了分潮的衰减,洪季潮波上溯至镇江处就基本消散。

5 摩阻作用及对余水位的控制

5.1 潮平均摩阻时空变化

根据Buschman等[4]推导的式(3)—式(5),计算洪枯季长江口主要站点径流作用、径潮相互作用、潮不对称作用对潮平均摩阻的贡献,结果见图4(其中S=Sr+Srt+St)。

图4 洪枯季径流作用、径潮相互作用、潮不对称作用对潮平均摩阻的贡献Fig.4 Contributions of river flow, river-tide interactions and tidal asymmetry to tidally averaged friction in wet season and dry season

从图4可以看出,长江口的潮平均摩阻以径流作用及径潮相互作用产生的摩阻为主,潮不对称作用的贡献微乎其微。究其原因,分潮簇流速振幅一般比日均流速小,加上算式内还包含分潮簇相位差的影响,使得St较小。

分析各摩阻分量的沿程变化规律,无论洪季还是枯季,Sr都从吴淞向上游逐渐增加,在南京取得最大值。产生这种现象的原因是:自上游向下传播过程中,由于和潮汐的相互作用以及河道断面的放宽,径流逐渐减弱,因此径流作用对潮平均摩阻的贡献逐渐减小。Srt从吴淞向上游大体上逐渐减小,在枯季沿程减小较慢,在洪季减小较快,这与图3中分潮簇洪枯季流速振幅沿程变化相似,表明径流增加会加剧潮波的衰减,从而显著降低径潮相互作用对潮平均摩阻的贡献。

对比洪枯季各站点Sr和Srt两项,Sr均>Srt,说明在长江高径流的背景下径流作用对潮平均摩阻的贡献始终更大。在江阴及以上站点,Sr占据主导,而在下游的徐六泾和吴淞站点,Srt对潮平均摩阻的贡献也较大,不可忽略。

5.2 摩阻分量对余水位的拟合

在一维动量守恒方程中,余水位梯度项基本与潮平均摩阻项相平衡,摩阻分量与余水位梯度间大致呈线性关系,因此通过摩阻分量的线性回归可以预测拟合余水位的变化。由5.1节可知St较小,因而选取潮平均摩阻中主要的Sr和Srt两项,建立二元线性回归模型,在各站点根据二者的日均值预测拟合对应的余水位,结果见图5。

图5 二元线性回归模型预测的余水位和数值模型模拟的余水位对比Fig.5 Comparison of residual water level predicted by binary linear regression model and simulated by numerical model

在线性回归方程中,p值反映检验的显著性水平,而决定系数R2反映自变量对于因变量变化的解释程度,本文中即余水位变化被摩阻分量变化定性解释的程度。

计算表明,所选取的5个站点二元线性回归模型的p值均<0.01,说明检验具有显著统计学差异,拟合方程可靠。图5(a)—图5(e),从上游南京至下游吴淞各站点二元线性回归各模型的R2分别为0.958、0.955、0.891、0.711、0.453,除最下游的吴淞站点外,均拟合良好,说明Sr和Srt的变化能较好地解释余水位的变化。除Sr和Srt的二元线性回归外,采用R2的大小作为依据,对比Sr、Srt、S各自的一元线性回归拟合结果,如表2所示。

表2 一元、二元线性回归拟合对比Table 2 Comparison between univariate and binary linear regression fittings

从表2可看出,摩阻分量之和S的一元线性回归即可较好地预测余水位的变化,体现出摩阻控制水位。由此可对图5中洪季拟合结果较枯季更好作出解释。因为洪季潮平均摩阻更大,动量方程中其他项的影响显得相对较小,因此余水位变化能被摩阻变化更好地预测。同理,由于吴淞站点位置最靠下,在洪季径流作用也较弱,因而拟合结果比其上游的站点要差一些。由于Sr始终>Srt,因此前者的一元线性拟合结果更优。在江阴以下站点,Srt对总摩阻的贡献变大,因此考虑了Srt变化的二元线性回归的拟合结果较一元线性回归有一定提升。

根据上述分析对图2所示长江口余水位时空变化进行解释。由于潮平均摩阻一直存在,平衡着余水位梯度,因此自吴淞向上游余水位一直增大。洪季径流量巨大,径流作用产生的摩阻占据主导且远大于枯季,使得洪季余水位梯度更大,余水位向上游增幅明显。枯季径流量变化不大,余水位的大小潮变化来源于潮平均摩阻的大小潮变化,归根到底是外海潮差变化所引起。

6 结 论

本文基于长江口二维水动力模型获得主要站点一年的潮位和流速数据,分析余水位时空变化特征;利用连续小波变换方法对流速进行分解,得到主要分潮簇流速的调和常数;通过摩阻公式计算径流作用、径潮相互作用和潮不对称作用对潮平均摩阻的贡献,并结合简单线性回归模型分析径潮动力相互作用下余水位的时空变化。得到的主要结论如下:

(1)在径潮动力相互作用下,长江口余水位自下游向上游逐渐抬升,表现出明显的洪枯季变化和大小潮变化;洪季余水位沿程增幅较大,枯季余水位沿程增幅较小,大潮时余水位比小潮时更高。

(2)长江口的潮平均摩阻以径流作用和径潮相互作用产生的摩阻为主,潮不对称项的贡献可忽略,且洪季的摩阻大于枯季;因径流量巨大,各站点径流作用对潮平均摩阻的贡献始终超过径潮相互作用。

(3)潮平均摩阻控制余水位,通过建立线性回归模型用摩阻的变化较好地解释说明了余水位的变化,所采用的回归模型包括一元与二元,决定系数R2多达到0.7以上;长江余水位的变化主要由径流作用控制,在下游站点考虑了径潮相互作用的二元回归模型较一元模型拟合效果有一定提升,径潮相互作用的影响不可忽视。