黃春華,林中源*,鄒華志,黃鵬飛,許 偉

(1.水利部珠江河口治理與保護重點實驗室,廣東 廣州 510611;2.珠江水利委員會珠江水利科學研究院,廣東 廣州 510611)

在氣候變化和人類活動影響下,世界范圍內河口咸潮上溯都有明顯加劇的趨勢[1],咸潮運動規(guī)律及發(fā)生機理研究已成為河口三角洲水動力過程、水資源利用和水環(huán)境保護的熱點研究方向[2]。珠江河口是粵港澳大灣區(qū)的核心區(qū)域,在國家發(fā)展大局中具有重要戰(zhàn)略地位,珠江河口周邊城市供水以河道型水源為主,且部分區(qū)域調蓄能力不足,枯季供水安全極易受河口咸潮威脅。尤其在磨刀門河口,咸潮災害頻發(fā),對澳門特別行政區(qū)、珠海和中山供水安全造成了嚴重威脅,地理位置見圖1。

圖1 磨刀門地理位置以及各站點位置

河口咸潮上溯受徑流、潮汐、風和平均海平面等多種驅動力影響[3]。近些年針對磨刀門河口咸潮上溯這一問題,諸多學者展開了大量研究,對咸潮上溯運動規(guī)律與發(fā)生機理有了較為深入的認識[4-10]。磨刀門咸潮上溯距離一般與徑流大小成反比,并基本呈冪函數(shù)關系[5];咸潮上溯峰值發(fā)生于小潮后的中潮期[4-7];海平面上升導致咸潮上溯加劇[8];人類活動導致的河床下切使潮動力增強,從而加劇咸潮上溯[9-10]。但磨刀門咸潮上溯對風的響應未有系統(tǒng)研究,一些學者[11-13]認為偏北風有利于磨刀門咸潮上溯,但劉雪峰等[14]發(fā)現(xiàn)偏東風導致咸潮上溯加劇。一般而言,風的作用可分為兩部分,一是跨岸風(垂直于岸線的風),二是沿岸風(平行于岸線的風)?？绨讹L直接作用于河口水體,通過離岸風和向岸風的應變作用改變水體的分層-混合狀態(tài),離岸風加強層化,向岸風加強混合,從而對咸潮上溯產(chǎn)生影響[15-17]。沿岸風通過改變口門外水位和口門邊界處鹽度,從而間接影響咸潮上溯[18-19]。對于磨刀門河口枯季而言,跨岸風(偏北風)和沿岸風(偏東風)均會對咸潮上溯產(chǎn)生一定影響,出現(xiàn)以上風對咸潮上溯的影響觀點不一致的原因為:以上研究關注的是事件性或短時間尺度內咸潮上溯動態(tài)特征,未能對磨刀門咸潮較長時間尺度的變化特征展開全面、系統(tǒng)研究。王藝霖等[20]探討了1998—2015年徑流量、平均海平面、平均潮差、最高潮位、最低潮位以及降雨量等因素與超標時數(shù)之間的相關性,但風的作用卻未被探討。

基于此,本文以磨刀門水道咸潮上溯為研究對象,在統(tǒng)計分析2004—2016年磨刀門水道7個站點日含氯度超標時間變化趨勢的基礎上,選取高要站日均流量、三灶站日最大潮差、石壁站日平均海平面、區(qū)域平均跨岸風和沿岸風等影響因素,利用經(jīng)驗正交函數(shù)分析(EOF)、相關分析及一元和多元線性回歸分析的方法,探究日超標時間與各影響因素的響應關系,量化各影響因素對日超標時間變化的貢獻,旨在識別跨岸風與沿岸風對磨刀門咸潮上溯的作用大小,并揭示磨刀門咸潮上溯主控因素,為區(qū)域調度及咸潮預報提供理論支撐。

1 研究資料

含氯度超標時間、流量、潮汐、平均海平面站點及風的取值范圍分布見圖1。2004—2016年磨刀門水道各監(jiān)測站點日含氯度超標時間數(shù)據(jù)來自中山水務局,包括大涌口、燈籠山、聯(lián)石灣、馬角、南鎮(zhèn)水廠、西河水廠、全祿水廠等7個站點,日超標時間指當日咸潮超過供水水源氯度上限250 mg/L的時段。每年數(shù)據(jù)從12月1日到翌年2月28日,但2010—2011年數(shù)據(jù)有缺失,僅有2010年12月1日到2011年1月29日的數(shù)據(jù)。2014—2016年高要站徑流量三灶潮汐數(shù)據(jù)來自水利部珠江水利委員會。每小時風的數(shù)據(jù)來自美國國家環(huán)境預測中心氣候預報系統(tǒng)再分析數(shù)據(jù)集,本文采用Lin等[21]研究磨刀門咸潮上溯的風場選取范圍,范圍包括磨刀門水道下游以及伶仃洋西部部分區(qū)域,并將風分解為跨岸風與沿岸風,沿岸風設為與磨刀門河口外的沿岸方向一致(東偏北22.5°),而跨岸風方向設為與沿岸風方向垂直,見圖1。2004—2016年香港石壁站的日平均海平面數(shù)據(jù)來自香港天文臺。

2 研究方法

由于磨刀門水道共有7個監(jiān)測站點,各站間的超標時間變化比較劇烈,用一個站點的數(shù)據(jù)不能代表河口的整體超標時間變化,因此對7個站點2004—2016年的日超標時間數(shù)據(jù)進行經(jīng)驗正交函數(shù)分析(EOF)。EOF分析類似于主成分分析(PCA),所分析的數(shù)據(jù)可看作是空間坐標(x,y)和時間t的函數(shù)s(x,y,t)。EOF分析將函數(shù)s(x,y,t)分解為空間坐標中一系列的正交函數(shù)fi(x,y),以及代表時間變化的一系列的時間函數(shù)gi(t)的累加。通常fi(x,y)稱為空間模態(tài),gi(t)稱為時間模態(tài),見式(1):

(1)

其中N為觀測站點的總數(shù)。通過求解由s(x,y,t)的協(xié)方差矩陣構造的特征值方程,可以得到fi(x,y)空間模態(tài)和其各自的系數(shù)gi(t)時間模態(tài),正交函數(shù)按協(xié)方差矩陣對應特征值的降序排列。因此,前幾個正交函數(shù)通?？山忉屗治鰯?shù)據(jù)大部分的變化特征。

3 結果與討論

3.1 磨刀門7個站點日超標時間年際變化

2004—2016年磨刀門河口日超標時間的時間序列見圖2,其中站點1—7分別代表了從水道下游到上游分布的大涌口、燈籠山、聯(lián)石灣、馬角、南鎮(zhèn)水廠、西河水廠和全祿水廠。由圖可知,磨刀門咸潮上溯年際變化較大。除2015—2016年外,其他年份全祿水廠站均出現(xiàn)咸潮上溯情況。通過EOF分析得到表征2004—2016年日超標時間的第一時間模態(tài)與第一空間模態(tài),見圖3。第一模態(tài)可解釋87%的超標時間變化,因此第一時間模態(tài)(PC1)能較大程度地反映2004—2016年磨刀門咸潮上溯的時間變化,并與圖2所示的咸潮上溯變化高度一致,下文將用PC1來表征磨刀門河口咸潮上溯的時間變化。

圖2 2004—2016年磨刀門水道枯季7個站點日超標時間的時間序列

a)第一時間模態(tài)

3.2 磨刀門咸潮上溯與各影響因素的關系

圖4展示了PC1隨徑流量、日最大潮差、平均海平面和風的日均及月均變化規(guī)律。值得注意的是,2007年2月(2006—2007年第三個紅星點),當月徑流(圖4b)、潮差(圖4c)、風速(圖4e)、平均海平面(圖4d)較上一個月平均值有所減小時,PC1也隨之減小,此時風向由東北風轉為東南風。類似的結果還出現(xiàn)在2008年1月及2009年2月。以上月均結果表明,徑流與咸潮上溯并不一直保持負相關關系,在徑流維持較低水平時,小幅度的增加徑流并不能有效抑制咸情,風的作用此時較為顯著。

a)2004—2016年表征日超標時間的PC1時間序列

由于外力一般是不穩(wěn)定的,河口在不斷變化的外力作用時一般處于非穩(wěn)定狀態(tài),并存在著一定的時間滯后[22-24]。本文利用相關分析判斷磨刀門咸潮上溯對各影響因素的滯后時間,當最大潮差、徑流量、沿岸風、跨岸風和平均海平面與不同滯后時間的PC1間相關系數(shù)達到最大時,便將此滯后時間作為該影響因素與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年。咸潮上溯與潮差的滯后時間取決于徑流量大小,同樣的結果出現(xiàn)在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與各影響因素滯后時間平均值

續(xù)圖5 每日最大潮差、徑流、平均海平面、沿岸風、跨岸風與不同滯后時間的PC1之間的相關系數(shù)

圖5 每日最大潮差、徑流、平均海平面、沿岸風、跨岸風與不同滯后時間的PC1之間的相關系數(shù)

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——平均海平面,各變量已根據(jù)上文中PC1與各變量的滯后時間進行調整;a、b、c、d、e、f——回歸系數(shù)。

為了便于各不同單位、量值變量間的比較,本文將以上變量均按其變量范圍標準化(或稱為無維化)。標準化方法可見式(3):

(3)

一元及多元線性回歸結果見圖6,所有線性回歸的顯著性水平均小于0.01。縱坐標表示決定系數(shù),通過對單個影響因素或多個影響因素與PC1進行線性擬合得到,反映咸潮上溯變化能通過回歸關系被徑流、潮差、平均海平面、沿岸風與跨岸風解釋的比例大小。結果表明,各影響因素與PC1線性回歸之間的決定系數(shù)(r2)逐年變化,說明每年枯季影響咸潮上溯變化的主導因素均有所不同。總體來看,徑流和潮汐依然是磨刀門咸潮上溯的主導因素,同時風和平均海平面的作用也不能忽略。每年各影響因素與PC1線性擬合后得到的決定系數(shù)見表2。除2004—2005、2015—2016年多元線性回歸分析的決定系數(shù)較低外,其他年份各影響因素均能解釋咸潮超標時間日變化的60%以上,各因素間的非線性作用是導致咸潮上溯日變化的年平均貢獻只有60%以上的原因[26]。2015—2016年徑流量遠超年平均徑流量,枯季大部分時間段內氯度日超標時間均為0 h,波動較小,所以導致與各影響因素之間的決定系數(shù)偏低。考慮所有影響因素,多元線性回歸分析結果的年平均決定系數(shù)為0.62,剔除2004—2005、2015—2016年特殊情況外,各影響因素可解釋咸潮超標時間日變化的68.3%。總體來看,2004—2016年各影響因素對咸潮上溯的貢獻最大的是潮差,其次是徑流,兩者與磨刀門超標時間之間的決定系數(shù)分別為0.30和0.28,風的作用略大于平均海平面,兩者對咸潮上溯的貢獻分別為12%和11%,其中能導致河口分層-混合變化的跨岸風作用略大于導致磨刀門口門鹽度、流速和水位變化的沿岸風作用。

表2 PC1與各影響因素一元及多元線性回歸分析決定系數(shù)結果

a)Q

多元線性回歸結果顯示,回歸系數(shù)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年各影響因素對磨刀門咸潮日超標時間變化的貢獻最大的是潮差,其次是徑流,兩者與氯度超標時間之間的決定系數(shù)分別為0.30和0.28,風的作用略大于平均海平面,兩者對咸潮上溯的貢獻分別為12%和11%,其中能導致河口分層和混合加劇的跨岸風作用略大于導致磨刀門口門鹽度、流速和水位變化的沿岸風作用。

b)徑流和潮差增大會減弱磨刀門咸潮上溯,平均海平面上升和利于下降流發(fā)育的風(偏東風)增強導致咸潮上溯加劇,而咸潮上溯對跨岸風的響應在不斷變化,離岸風與向岸風均會造成咸潮上溯加劇。

c)通過統(tǒng)計2004—2016年磨刀門水道咸潮日超標時間與各影響因素的滯后時間得出,日超標時間滯后于潮差3.2 d,滯后于徑流1.8 d,滯后于平均海平面1.0 d,滯后于沿岸風、跨岸風1.2 d。意味著磨刀門咸潮峰值出現(xiàn)在小潮潮差最小后的3.2 d,而高要流量增大,1.8 d后才能抑制咸情。本文可為區(qū)域調度及咸潮預報提供理論支撐。