磨刀门枯季咸潮上溯归因分析与响应规律研究

黄春华,林中源*,邹华志,黄鹏飞,许 伟

(1.水利部珠江河口治理与保护重点实验室,广东 广州 510611;2.珠江水利委员会珠江水利科学研究院,广东 广州 510611)

在气候变化和人类活动影响下,世界范围内河口咸潮上溯都有明显加剧的趋势[1],咸潮运动规律及发生机理研究已成为河口三角洲水动力过程、水资源利用和水环境保护的热点研究方向[2]。珠江河口是粤港澳大湾区的核心区域,在国家发展大局中具有重要战略地位,珠江河口周边城市供水以河道型水源为主,且部分区域调蓄能力不足,枯季供水安全极易受河口咸潮威胁。尤其在磨刀门河口,咸潮灾害频发,对澳门特别行政区、珠海和中山供水安全造成了严重威胁,地理位置见图1。

图1 磨刀门地理位置以及各站点位置

河口咸潮上溯受径流、潮汐、风和平均海平面等多种驱动力影响[3]。近些年针对磨刀门河口咸潮上溯这一问题,诸多学者展开了大量研究,对咸潮上溯运动规律与发生机理有了较为深入的认识[4-10]。磨刀门咸潮上溯距离一般与径流大小成反比,并基本呈幂函数关系[5];咸潮上溯峰值发生于小潮后的中潮期[4-7];海平面上升导致咸潮上溯加剧[8];人类活动导致的河床下切使潮动力增强,从而加剧咸潮上溯[9-10]。但磨刀门咸潮上溯对风的响应未有系统研究,一些学者[11-13]认为偏北风有利于磨刀门咸潮上溯,但刘雪峰等[14]发现偏东风导致咸潮上溯加剧。一般而言,风的作用可分为两部分,一是跨岸风(垂直于岸线的风),二是沿岸风(平行于岸线的风)。跨岸风直接作用于河口水体,通过离岸风和向岸风的应变作用改变水体的分层-混合状态,离岸风加强层化,向岸风加强混合,从而对咸潮上溯产生影响[15-17]。沿岸风通过改变口门外水位和口门边界处盐度,从而间接影响咸潮上溯[18-19]。对于磨刀门河口枯季而言,跨岸风(偏北风)和沿岸风(偏东风)均会对咸潮上溯产生一定影响,出现以上风对咸潮上溯的影响观点不一致的原因为:以上研究关注的是事件性或短时间尺度内咸潮上溯动态特征,未能对磨刀门咸潮较长时间尺度的变化特征展开全面、系统研究。王艺霖等[20]探讨了1998—2015年径流量、平均海平面、平均潮差、最高潮位、最低潮位以及降雨量等因素与超标时数之间的相关性,但风的作用却未被探讨。

基于此,本文以磨刀门水道咸潮上溯为研究对象,在统计分析2004—2016年磨刀门水道7个站点日含氯度超标时间变化趋势的基础上,选取高要站日均流量、三灶站日最大潮差、石壁站日平均海平面、区域平均跨岸风和沿岸风等影响因素,利用经验正交函数分析(EOF)、相关分析及一元和多元线性回归分析的方法,探究日超标时间与各影响因素的响应关系,量化各影响因素对日超标时间变化的贡献,旨在识别跨岸风与沿岸风对磨刀门咸潮上溯的作用大小,并揭示磨刀门咸潮上溯主控因素,为区域调度及咸潮预报提供理论支撑。

1 研究资料

含氯度超标时间、流量、潮汐、平均海平面站点及风的取值范围分布见图1。2004—2016年磨刀门水道各监测站点日含氯度超标时间数据来自中山水务局,包括大涌口、灯笼山、联石湾、马角、南镇水厂、西河水厂、全禄水厂等7个站点,日超标时间指当日咸潮超过供水水源氯度上限250 mg/L的时段。每年数据从12月1日到翌年2月28日,但2010—2011年数据有缺失,仅有2010年12月1日到2011年1月29日的数据。2014—2016年高要站径流量三灶潮汐数据来自水利部珠江水利委员会。每小时风的数据来自美国国家环境预测中心气候预报系统再分析数据集,本文采用Lin等[21]研究磨刀门咸潮上溯的风场选取范围,范围包括磨刀门水道下游以及伶仃洋西部部分区域,并将风分解为跨岸风与沿岸风,沿岸风设为与磨刀门河口外的沿岸方向一致(东偏北22.5°),而跨岸风方向设为与沿岸风方向垂直,见图1。2004—2016年香港石壁站的日平均海平面数据来自香港天文台。

2 研究方法

由于磨刀门水道共有7个监测站点,各站间的超标时间变化比较剧烈,用一个站点的数据不能代表河口的整体超标时间变化,因此对7个站点2004—2016年的日超标时间数据进行经验正交函数分析(EOF)。EOF分析类似于主成分分析(PCA),所分析的数据可看作是空间坐标(x,y)和时间t的函数s(x,y,t)。EOF分析将函数s(x,y,t)分解为空间坐标中一系列的正交函数fi(x,y),以及代表时间变化的一系列的时间函数gi(t)的累加。通常fi(x,y)称为空间模态,gi(t)称为时间模态,见式(1):

(1)

其中N为观测站点的总数。通过求解由s(x,y,t)的协方差矩阵构造的特征值方程,可以得到fi(x,y)空间模态和其各自的系数gi(t)时间模态,正交函数按协方差矩阵对应特征值的降序排列。因此,前几个正交函数通常可解释所分析数据大部分的变化特征。

3 结果与讨论

3.1 磨刀门7个站点日超标时间年际变化

2004—2016年磨刀门河口日超标时间的时间序列见图2,其中站点1—7分别代表了从水道下游到上游分布的大涌口、灯笼山、联石湾、马角、南镇水厂、西河水厂和全禄水厂。由图可知,磨刀门咸潮上溯年际变化较大。除2015—2016年外,其他年份全禄水厂站均出现咸潮上溯情况。通过EOF分析得到表征2004—2016年日超标时间的第一时间模态与第一空间模态,见图3。第一模态可解释87%的超标时间变化,因此第一时间模态(PC1)能较大程度地反映2004—2016年磨刀门咸潮上溯的时间变化,并与图2所示的咸潮上溯变化高度一致,下文将用PC1来表征磨刀门河口咸潮上溯的时间变化。

图2 2004—2016年磨刀门水道枯季7个站点日超标时间的时间序列

a)第一时间模态

3.2 磨刀门咸潮上溯与各影响因素的关系

图4展示了PC1随径流量、日最大潮差、平均海平面和风的日均及月均变化规律。值得注意的是,2007年2月(2006—2007年第三个红星点),当月径流(图4b)、潮差(图4c)、风速(图4e)、平均海平面(图4d)较上一个月平均值有所减小时,PC1也随之减小,此时风向由东北风转为东南风。类似的结果还出现在2008年1月及2009年2月。以上月均结果表明,径流与咸潮上溯并不一直保持负相关关系,在径流维持较低水平时,小幅度的增加径流并不能有效抑制咸情,风的作用此时较为显著。

a)2004—2016年表征日超标时间的PC1时间序列

由于外力一般是不稳定的,河口在不断变化的外力作用时一般处于非稳定状态,并存在着一定的时间滞后[22-24]。本文利用相关分析判断磨刀门咸潮上溯对各影响因素的滞后时间,当最大潮差、径流量、沿岸风、跨岸风和平均海平面与不同滞后时间的PC1间相关系数达到最大时,便将此滞后时间作为该影响因素与PC1的最终滞后时间。图5所示,不同年份PC1与不同影响因素之间的滞后时间存在年际变化。在低径流量(<2 000 m3/s)时,PC1滞后于日最大潮差4 d,如2004—2005、2005—2006、2007—2008、2011—2012年;而在高径流量(>2 200 m3/s)时,PC1滞后于日最大潮差3 d,如2006—2007、2008—2009、2012—2013、2013—2014、2014—2015年。咸潮上溯与潮差的滞后时间取决于径流量大小,同样的结果出现在PC1与平均海平面、跨岸风、沿岸风的滞后时间上。当径流量较大时,咸潮上溯对该3种影响因素的滞后时间较短;径流量较小时,咸潮上溯与该3种影响因素的滞后时间较长。2004—2016年可代表磨刀门水道整体超标时间的PC1与各影响因素不同滞后时间平均值结果见表1。结果表明:磨刀门水道咸潮上溯滞后于潮差3.2 d,滞后于径流1.8 d,滞后于平均海平面1.0 d,滞后于沿岸风、跨岸风1.2 d。咸潮上溯与潮差、径流的滞后时间,与Gong等[5]及Liu等[25]基本一致。

表1 2004—2016年PC1与各影响因素滞后时间平均值

续图5 每日最大潮差、径流、平均海平面、沿岸风、跨岸风与不同滞后时间的PC1之间的相关系数

图5 每日最大潮差、径流、平均海平面、沿岸风、跨岸风与不同滞后时间的PC1之间的相关系数

3.3 各影响因素对咸潮上溯的影响及贡献

为了探讨2004—2016年咸潮上溯的主控因素及其他各影响因素对咸潮上溯的贡献大小,本文利用一元和多元线性回归进行分析,并将上文中每年PC1与各影响因素之间的滞后时间应用到回归分析当中。

回归分析中考虑的影响因素为径流、潮差、平均海平面、沿岸风以及跨岸风,见式(2):

S=a+b×Q+c×T+d×U+e×V+f×M

(2)

式中S——PC1;Q——日径流量;T——日最大潮差;U——日沿岸风风应力;V——日跨岸风风应力;M——平均海平面,各变量已根据上文中PC1与各变量的滞后时间进行调整;a、b、c、d、e、f——回归系数。

为了便于各不同单位、量值变量间的比较,本文将以上变量均按其变量范围标准化(或称为无维化)。标准化方法可见式(3):

(3)

一元及多元线性回归结果见图6,所有线性回归的显著性水平均小于0.01。纵坐标表示决定系数,通过对单个影响因素或多个影响因素与PC1进行线性拟合得到,反映咸潮上溯变化能通过回归关系被径流、潮差、平均海平面、沿岸风与跨岸风解释的比例大小。结果表明,各影响因素与PC1线性回归之间的决定系数(r2)逐年变化,说明每年枯季影响咸潮上溯变化的主导因素均有所不同。总体来看,径流和潮汐依然是磨刀门咸潮上溯的主导因素,同时风和平均海平面的作用也不能忽略。每年各影响因素与PC1线性拟合后得到的决定系数见表2。除2004—2005、2015—2016年多元线性回归分析的决定系数较低外,其他年份各影响因素均能解释咸潮超标时间日变化的60%以上,各因素间的非线性作用是导致咸潮上溯日变化的年平均贡献只有60%以上的原因[26]。2015—2016年径流量远超年平均径流量,枯季大部分时间段内氯度日超标时间均为0 h,波动较小,所以导致与各影响因素之间的决定系数偏低。考虑所有影响因素,多元线性回归分析结果的年平均决定系数为0.62,剔除2004—2005、2015—2016年特殊情况外,各影响因素可解释咸潮超标时间日变化的68.3%。总体来看,2004—2016年各影响因素对咸潮上溯的贡献最大的是潮差,其次是径流,两者与磨刀门超标时间之间的决定系数分别为0.30和0.28,风的作用略大于平均海平面,两者对咸潮上溯的贡献分别为12%和11%,其中能导致河口分层-混合变化的跨岸风作用略大于导致磨刀门口门盐度、流速和水位变化的沿岸风作用。

表2 PC1与各影响因素一元及多元线性回归分析决定系数结果

a)Q

多元线性回归结果显示,回归系数b、c、d为负,f为正,而e在2007—2008、2010—2011、2012—2013、2013—2014年为正,其他年份为负,与图5中PC1与各影响因素相关关系保持一致。以上结果表明,径流和潮差增大会减弱咸潮上溯,平均海平面上升和沿岸风(偏东风)增强导致咸潮上溯加剧,而咸潮上溯对跨岸风的响应在不断变化,离岸风与向岸风均会造成咸潮上溯加剧,与以往研究结果不一致[11-13]。本文推测,与2007—2008、2010—2011、2012—2013、2013—2014年径流均有一个较大的波动(图4b)有关,咸潮上溯主要随着径流量波动而变化,从而导致跨岸风的变化与咸潮上溯变化呈正相关;或与当时河口层化较弱有关,向岸风通过平流作用将外海高盐水输送到磨刀门上游,具体原因值得深入探讨。

4 结论

利用EOF分析、相关性分析以及一元和多元线性回归分析等分析方法,研究了2004—2016年磨刀门河口枯季咸潮上溯变化特征及主控因素,结论如下。

a)2004—2016年各影响因素对磨刀门咸潮日超标时间变化的贡献最大的是潮差,其次是径流,两者与氯度超标时间之间的决定系数分别为0.30和0.28,风的作用略大于平均海平面,两者对咸潮上溯的贡献分别为12%和11%,其中能导致河口分层和混合加剧的跨岸风作用略大于导致磨刀门口门盐度、流速和水位变化的沿岸风作用。

b)径流和潮差增大会减弱磨刀门咸潮上溯,平均海平面上升和利于下降流发育的风(偏东风)增强导致咸潮上溯加剧,而咸潮上溯对跨岸风的响应在不断变化,离岸风与向岸风均会造成咸潮上溯加剧。

c)通过统计2004—2016年磨刀门水道咸潮日超标时间与各影响因素的滞后时间得出,日超标时间滞后于潮差3.2 d,滞后于径流1.8 d,滞后于平均海平面1.0 d,滞后于沿岸风、跨岸风1.2 d。意味着磨刀门咸潮峰值出现在小潮潮差最小后的3.2 d,而高要流量增大,1.8 d后才能抑制咸情。本文可为区域调度及咸潮预报提供理论支撑。