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

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

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

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

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

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

1 研究資料

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

2 研究方法

由于磨刀門水道共有7個監(jiān)測站點,各站間的超標(biāo)時間變化比較劇烈,用一個站點的數(shù)據(jù)不能代表河口的整體超標(biāo)時間變化,因此對7個站點2004—2016年的日超標(biāo)時間數(shù)據(jù)進(jìn)行經(jīng)驗正交函數(shù)分析(EOF)。EOF分析類似于主成分分析(PCA),所分析的數(shù)據(jù)可看作是空間坐標(biāo)(x,y)和時間t的函數(shù)s(x,y,t)。EOF分析將函數(shù)s(x,y,t)分解為空間坐標(biāo)中一系列的正交函數(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é)方差矩陣構(gòu)造的特征值方程,可以得到fi(x,y)空間模態(tài)和其各自的系數(shù)gi(t)時間模態(tài),正交函數(shù)按協(xié)方差矩陣對應(yīng)特征值的降序排列。因此,前幾個正交函數(shù)通?？山忉屗治鰯?shù)據(jù)大部分的變化特征。

3 結(jié)果與討論

3.1 磨刀門7個站點日超標(biāo)時間年際變化

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

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

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

3.2 磨刀門咸潮上溯與各影響因素的關(guān)系

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

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

由于外力一般是不穩(wěn)定的,河口在不斷變化的外力作用時一般處于非穩(wěn)定狀態(tài),并存在著一定的時間滯后[22-24]。本文利用相關(guān)分析判斷磨刀門咸潮上溯對各影響因素的滯后時間,當(dāng)最大潮差、徑流量、沿岸風(fēng)、跨岸風(fēng)和平均海平面與不同滯后時間的PC1間相關(guān)系數(shù)達(dá)到最大時,便將此滯后時間作為該影響因素與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年。咸潮上溯與潮差的滯后時間取決于徑流量大小,同樣的結(jié)果出現(xiàn)在PC1與平均海平面、跨岸風(fēng)、沿岸風(fēng)的滯后時間上。當(dāng)徑流量較大時,咸潮上溯對該3種影響因素的滯后時間較短;徑流量較小時,咸潮上溯與該3種影響因素的滯后時間較長。2004—2016年可代表磨刀門水道整體超標(biāo)時間的PC1與各影響因素不同滯后時間平均值結(jié)果見表1。結(jié)果表明:磨刀門水道咸潮上溯滯后于潮差3.2 d,滯后于徑流1.8 d,滯后于平均海平面1.0 d,滯后于沿岸風(fēng)、跨岸風(fēng)1.2 d。咸潮上溯與潮差、徑流的滯后時間,與Gong等[5]及Liu等[25]基本一致。

表1 2004—2016年P(guān)C1與各影響因素滯后時間平均值

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

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

3.3 各影響因素對咸潮上溯的影響及貢獻(xiàn)

為了探討2004—2016年咸潮上溯的主控因素及其他各影響因素對咸潮上溯的貢獻(xiàn)大小,本文利用一元和多元線性回歸進(jìn)行分析,并將上文中每年P(guān)C1與各影響因素之間的滯后時間應(yīng)用到回歸分析當(dāng)中。

回歸分析中考慮的影響因素為徑流、潮差、平均海平面、沿岸風(fēng)以及跨岸風(fēng),見式(2):

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

(2)

式中S——PC1;Q——日徑流量;T——日最大潮差;U——日沿岸風(fēng)風(fēng)應(yīng)力;V——日跨岸風(fēng)風(fēng)應(yīng)力;M——平均海平面,各變量已根據(jù)上文中PC1與各變量的滯后時間進(jìn)行調(diào)整;a、b、c、d、e、f——回歸系數(shù)。

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

(3)

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

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

a)Q

多元線性回歸結(jié)果顯示,回歸系數(shù)b、c、d為負(fù),f為正,而e在2007—2008、2010—2011、2012—2013、2013—2014年為正,其他年份為負(fù),與圖5中PC1與各影響因素相關(guān)關(guān)系保持一致。以上結(jié)果表明,徑流和潮差增大會減弱咸潮上溯,平均海平面上升和沿岸風(fēng)(偏東風(fēng))增強導(dǎo)致咸潮上溯加劇,而咸潮上溯對跨岸風(fēng)的響應(yīng)在不斷變化,離岸風(fēng)與向岸風(fēng)均會造成咸潮上溯加劇,與以往研究結(jié)果不一致[11-13]。本文推測,與2007—2008、2010—2011、2012—2013、2013—2014年徑流均有一個較大的波動(圖4b)有關(guān),咸潮上溯主要隨著徑流量波動而變化,從而導(dǎo)致跨岸風(fēng)的變化與咸潮上溯變化呈正相關(guān);或與當(dāng)時河口層化較弱有關(guān),向岸風(fēng)通過平流作用將外海高鹽水輸送到磨刀門上游,具體原因值得深入探討。

4 結(jié)論

利用EOF分析、相關(guān)性分析以及一元和多元線性回歸分析等分析方法,研究了2004—2016年磨刀門河口枯季咸潮上溯變化特征及主控因素,結(jié)論如下。

a)2004—2016年各影響因素對磨刀門咸潮日超標(biāo)時間變化的貢獻(xiàn)最大的是潮差,其次是徑流,兩者與氯度超標(biāo)時間之間的決定系數(shù)分別為0.30和0.28,風(fēng)的作用略大于平均海平面,兩者對咸潮上溯的貢獻(xiàn)分別為12%和11%,其中能導(dǎo)致河口分層和混合加劇的跨岸風(fēng)作用略大于導(dǎo)致磨刀門口門鹽度、流速和水位變化的沿岸風(fēng)作用。

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

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