珠江磨刀門河口水位分布演變特征及其對(duì)人類活動(dòng)的響應(yīng)*

馬玉婷, 蔡華陽, 楊昊, 劉鋒, 陳歐, 謝榮耀,歐素英, 楊清書

1. 中山大學(xué), 海洋工程與技術(shù)學(xué)院, 河口海岸研究所, 廣東 廣州 510275;

2. 河口水利技術(shù)國家地方聯(lián)合工程實(shí)驗(yàn)室, 廣東 廣州 510275;

3. 廣東省海岸與島礁工程技術(shù)研究中心, 廣東 廣州 510275;

4. 南方海洋科學(xué)與工程廣東省實(shí)驗(yàn)室(珠海), 廣東 珠海 519000

河口是海陸相互作用的過渡區(qū)域, 受徑流和潮流的相互作用。余水位(即潮平均水位)是河口徑流、潮汐、地形等多因素耦合的結(jié)果, 其空間形態(tài)變化及演變機(jī)制是探討河口陸海相互作用及潮波衰減的重要基礎(chǔ)及前沿問題, 對(duì)探討河口區(qū)水資源高效開發(fā)利用具有重要科學(xué)意義。河口區(qū)潮平均狀態(tài)下,其水位坡降(水位空間變化率)近似與底摩擦項(xiàng)(或底切應(yīng)力項(xiàng))相平衡(Buschman et al, 2009; Cai et al,2015; Hoitink et al, 2016), 因此平均態(tài)下三角洲河道及河口區(qū)余水面坡降梯度(即余水位曲率)與河床底切應(yīng)力梯度相平衡, 可指示底切應(yīng)力梯度的時(shí)空演變, 并在一定程度上指示水道泥沙凈輸運(yùn)的方向和地形沖淤趨勢(shì)。研究河口余水位的空間變化尤其是余水位曲率的時(shí)空變化, 可從水動(dòng)力角度揭示自然因素和人類活動(dòng)共同影響下河床沖淤演變趨勢(shì),為防洪納潮、水沙調(diào)控和航道整治開發(fā)等提供重要理論支撐。

密度梯度、潮汐不對(duì)稱和回水效應(yīng)的共同作用使得河口的平均水位向陸地方向抬升, 余水位及其坡降隨流量增大而增大(Cai et al, 2013, 2014)。但余水位坡降隨流量的增大不是單調(diào)遞增, 而是存在閾值效應(yīng)(張先毅 等, 2019)。低流量時(shí)河口過渡區(qū)呈凹形水面輪廓, 形成空間減速和沉積區(qū), 而高流量時(shí)形成水流加速和侵蝕區(qū)(Lamb et al, 2012)。近年來,強(qiáng)人類活動(dòng)對(duì)河口的干預(yù)日趨嚴(yán)重, 超過自然恢復(fù)能力, 已成為河口三角洲演變過程中的第三驅(qū)動(dòng)力(陳吉余 等, 2008)。聯(lián)圍筑閘、灘涂圍墾、河道采沙等強(qiáng)人類活動(dòng)已經(jīng)對(duì)三角洲河床演變產(chǎn)生重要影響。作為珠江主要泄洪通道的磨刀門河口, 近年來在挖沙等人類活動(dòng)影響下上游河道河床平均加深約0.7m, 中下游河道河床加深約1.4m, 河道逐漸向窄深化方向發(fā)展(錢挹清, 2004; 鄭國棟, 2005; 劉鋒等, 2011)。河床下切導(dǎo)致河口上段水位顯著下降、潮汐動(dòng)力增強(qiáng)、咸潮入侵加劇等(Cai et al, 2019; 洪鵬鋒 等, 2019)。

水位作為河口三角洲地區(qū)觀測(cè)站位最多、時(shí)間序列最長且最連續(xù)的水文要素, 其在自然和人類活動(dòng)影響下的時(shí)間演變特征及機(jī)制都得到較深入的研究。例如楊昊等(2020)、張先毅等(2020)對(duì)磨刀門河口平均水位變化及機(jī)制進(jìn)行研究, 指出余水位在時(shí)間演變上存在著明顯的季節(jié)變化和人類活動(dòng)影響下的階段性異變特征, 并采用邊界驅(qū)動(dòng)的回歸模型對(duì)挖沙等強(qiáng)人類活動(dòng)的影響進(jìn)行定量辨識(shí), 指出高強(qiáng)度采沙導(dǎo)致的河床下切使河口上段水位與水面坡降顯著減小, 特別是夏季洪水流量時(shí)水面坡降減小最大。但針對(duì)余水位的空間分布及其演變特征、水位曲率變化對(duì)強(qiáng)人類活動(dòng)的響應(yīng)以及其與地形沖淤演變的關(guān)聯(lián)性等問題仍有待進(jìn)一步深入研究。本文聚焦磨刀門河口余水位空間分布及其演變特征, 基于長序列水位觀測(cè)數(shù)據(jù), 在探究余水位空間分布的基礎(chǔ)上, 采用雙變量線性回歸模型定量辨識(shí)人類活動(dòng)對(duì)余水位坡降和曲率的影響, 初步從水動(dòng)力演變角度探討三角洲河道沖淤演變趨勢(shì)。

1 研究區(qū)域概況

珠江三角洲上承三江(西江、北江和東江), 下連八大口門(虎門、蕉門、洪奇門、橫門、磨刀門、雞啼門、虎跳門、崖門), 河網(wǎng)縱橫, 洪潮疊加, 是世界上最為復(fù)雜的河網(wǎng)三角洲體系之一。西江干流進(jìn)入珠江三角洲后于思賢 滘出現(xiàn)分流分沙, 其中思賢滘 西 滘口至百頃頭河段稱為西江干流水道(約80km),百頃頭至口門河段稱為磨刀門水道(約64km)。本文研究區(qū)域包括西江干流水道和磨刀門水道(以下簡(jiǎn)稱為磨刀門河口), 即從上游的馬口站延伸至下游口門附近的三灶站, 全程約144km (如圖1 所示)。磨刀門河口為西北江河網(wǎng)的主要泄洪通道, 其多年輸沙輸水量居八大口門之首, 屬于典型的河優(yōu)型河口(呂海濱 等, 2006)。據(jù)馬口水文站1965—2016 年月均流量資料統(tǒng)計(jì), 其多年洪季平均流量為10688m3·s-1, 多年枯季平均流量為3583m3·s-1, 洪枯季差異顯著。上游馬口水文站多年月均水位為1.43m, 月平均潮差為0.34m, 下游三灶潮位站多年月均水位為-0.10m, 月平均潮差為1.09m, 口門至上游沿程水位抬升, 潮差減小。20 世紀(jì)80 年代中期開始, 珠江三角洲河網(wǎng)區(qū)出現(xiàn)大規(guī)模的河道挖沙現(xiàn)象。據(jù)前人統(tǒng)計(jì), 西江干流水道1991—1997 年年均采沙量約為2.21×106m3, 采沙總深0.71m; 竹洲頭至燈籠山1991—2000 年年均采沙量約為3.16×106m3,采沙總深1.06m; 年均采沙量約為年均來沙量的5倍, 引起河床高程不均勻下切0.7～2.3m, 河道逐漸向窄深化方向發(fā)展(錢挹清, 2004; 劉鋒 等, 2011)。20 世紀(jì)60 年代中期起, 磨刀門口門地區(qū)開始大規(guī)?？陂T整治和圍墾工程, 1966—1988 年建設(shè)白藤堵海大堤, 1984—1995 年建設(shè)橫洲水道東西導(dǎo)堤和洪灣水道南北導(dǎo)堤, 口門區(qū)大片水域圍墾成陸, 其水域面積減小約80%, 入海口門向外海延伸約16km,口門水位逐漸抬升(黎開志 等, 2005; 韓志遠(yuǎn) 等,2010; 胡煌昊 等, 2016)。磨刀門主干道疏浚工程于1999 年啟動(dòng), 整治工程河段過洪斷面面積加大, 附近水位略有下降(袁建國 等, 2009)。受地形邊界(挖沙引起地形下切和口門圍墾引起河道變窄以及疏浚引起斷面面積加大等)和動(dòng)力邊界(海平面和上游流量變化等)的雙重影響, 磨刀門河口水位呈現(xiàn)其獨(dú)特的分布及階段性變化。

圖1 研究區(qū)域圖及測(cè)站位置該圖基于國家測(cè)繪地理信息局標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站下載的審圖號(hào)為GS(2019)4342 標(biāo)準(zhǔn)地圖制作Fig. 1 Map of the study area and the stations used in this study

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)

本文收集馬口、甘竹、江門、竹銀、燈籠山和三灶6 個(gè)測(cè)站1965—2016 年月平均水位數(shù)據(jù)以及相應(yīng)時(shí)段馬口站的月均流量數(shù)據(jù), 站位位置如圖1 所示。所用數(shù)據(jù)來源于廣東省水文局, 原始數(shù)據(jù)高程基面為凍結(jié)基面, 已統(tǒng)一校正至珠江基面高程, 基面數(shù)據(jù)缺失年份使用已知相鄰年份的基面進(jìn)行轉(zhuǎn)換,可能導(dǎo)致水位有0～0.07m 的偏差, 相對(duì)于河段長度屬于小量。

2.2 水面線沿程變化的數(shù)理統(tǒng)計(jì)

一維動(dòng)量守恒方程中, 水位坡降項(xiàng)是河道水流運(yùn)動(dòng)的主要驅(qū)動(dòng)力, 它是水位沿空間某個(gè)方向的變化率, 比如在一維河道中, 其水面坡降項(xiàng)是水位沿河道延伸方向, 即x方向(河道延伸方向)的一階導(dǎo)數(shù)。但實(shí)際上由于河道地形、徑潮動(dòng)力的復(fù)雜性, 無法得到河口三角洲地區(qū)水位坡降(S)的理論解, 因此實(shí)際應(yīng)用中主要采用其離散解來確定, 即:

由于在西江磨刀門河口144km 的范圍內(nèi), 實(shí)測(cè)站點(diǎn)有限, 為更直觀地顯示水位空間分布, 假設(shè)三灶—馬口之間的水位沿程變化可由此6 個(gè)站點(diǎn)(圖2b*號(hào))控制, 以三灶為起點(diǎn), 設(shè)其x=0, 河道余水位(Z)沿河道距離x的變化可分別使用分段三次Hermite 插值、四階多項(xiàng)式擬合、三次樣條插值進(jìn)行計(jì)算, 結(jié)果如圖2b 所示(圖中只為展示3 種方法的差異性, 其水位值非真實(shí)水位)。由圖2b 可見, 分段三次Hermite 插值方法可避免過沖, 且可準(zhǔn)確地連接平臺(tái)區(qū), 保證原始數(shù)據(jù)的準(zhǔn)確性, 更能反映出相鄰站位之間平均信息。因此, 本文主要采用分段三次Hermite 插值, 基于6 個(gè)站位的距離x與月平均余水位Z, 按空間間隔1km 對(duì)整個(gè)河口進(jìn)行插值, 得出沿程 144 個(gè)斷面的余水位Z, 并采用公式(1)和公式(2)分別計(jì)算沿程余水位坡降和余水位曲率, 用以定量分析月平均水位的空間分布及其變化規(guī)律。

圖2 余水位曲率示意圖(a)及水位沿程分布擬合方法對(duì)比(b)Fig. 2 Sketch of the residual water level curvature (a) and comparison of fitting methods for water level profiles along the estuary (b)

2.3 定量辨識(shí)強(qiáng)人類活動(dòng)影響的雙變量線性回歸模型

在三角洲地區(qū), 強(qiáng)人類活動(dòng)(如圍墾、大規(guī)模挖沙、河道疏浚等)導(dǎo)致三角洲地形發(fā)生顯著變化, 進(jìn)而引起徑潮動(dòng)力結(jié)構(gòu)的改變。Cai 等(2019)、Yang等(2020)采用上下動(dòng)力邊界驅(qū)動(dòng)的線性回歸模型,定量辨識(shí)河道下切等對(duì)三角洲徑潮動(dòng)力的影響。本文在此基礎(chǔ)上, 同樣采用基于數(shù)據(jù)驅(qū)動(dòng)的雙變量線性回歸模型, 以探討水位空間結(jié)構(gòu)的影響因素, 通過輸入流量(Q)和下游邊界水位(Zdown), 重構(gòu)出河口沿程各站點(diǎn)的平均水位(Zsim), 回歸模型的效果用均方根誤差(RMSE)和相關(guān)系數(shù)(R2)來評(píng)價(jià):

式中:α1、α2和β分別為流量、下游水位邊界的待定系數(shù)和背景值待定系數(shù)。n為輸入數(shù)據(jù)的長度,Zobs和Zsim分別表示實(shí)測(cè)水位和計(jì)算水位值, 上劃線表示平均值。

人類活動(dòng)對(duì)河口動(dòng)力格局的影響是一個(gè)不斷累積的過程, 累積到一定程度動(dòng)力格局會(huì)發(fā)生異變,根據(jù)Yang 等(2020)、張先毅等(2020)的研究, 磨刀門河口動(dòng)力結(jié)構(gòu)對(duì)人類活動(dòng)的響應(yīng)主要集中在1990—2000 年期間, 因此將1965 年至1990 年(共26年)定義為緩變階段(即強(qiáng)人類活動(dòng)前), 2001 年至2016 年(共16 年)定義為強(qiáng)人類活動(dòng)影響后的調(diào)整階段?；诟髡军c(diǎn)強(qiáng)人類活動(dòng)前的動(dòng)力邊界驅(qū)動(dòng)回歸模型及其重構(gòu)的月平均水位, 計(jì)算出強(qiáng)人類活動(dòng)前的沿程余水位坡降和余水位曲率; 使用緩變階段的待定系數(shù), 假設(shè)地形邊界條件不變, 采用調(diào)整階段的流量及水位邊界驅(qū)動(dòng), 預(yù)報(bào)并計(jì)算調(diào)整階段的余水位坡降和曲率。本文聚焦強(qiáng)人類活動(dòng)對(duì)余水位曲率變化的影響, 強(qiáng)人類活動(dòng)前后曲率總變化( ΔTOT)可分解為地形邊界影響( ΔGEO, 即挖沙等強(qiáng)人類活動(dòng)改變河道地形邊界后的影響)和動(dòng)力邊界影響( ΔBOU, 即上游徑流量和外海潮汐動(dòng)力改變后的影響)(Cai et al, 2019; 楊昊 等, 2020), 通過計(jì)算 ΔGEO來定量辨識(shí)強(qiáng)人類活動(dòng)對(duì)曲率的影響程度:

式中:Cobs和Csim分別為實(shí)測(cè)曲率值和計(jì)算曲率值,下標(biāo)pre 和post 分別表示緩變階段(強(qiáng)人類活動(dòng)前階段)和調(diào)整階段(強(qiáng)人類活動(dòng)后階段)

2.4 余水位曲率與底床沖淤趨勢(shì)的關(guān)聯(lián)性

式中:α為待定系數(shù)。

基于一維泥沙守恒方程, 河道底床沖淤變化模型可描述為:

3 水位分布演變特征

從實(shí)測(cè)水位數(shù)據(jù)出發(fā), 對(duì)沿程6 個(gè)站點(diǎn)月平均余水位進(jìn)行三次Hermite 插值(1998 至2000 年江門站無水位數(shù)據(jù), 故不進(jìn)行插值), 計(jì)算出沿程余水位坡降與余水位曲率的數(shù)據(jù)序列, 結(jié)果如圖3 所示。由圖3a 可見, 余水位從口門至上游沿程抬升, 上游洪枯季變化較下游顯著。由圖3b 可見, 在江門至甘竹段余水位坡降存在負(fù)向波動(dòng), 且這種波動(dòng)多在枯季出現(xiàn), 而在強(qiáng)人類活動(dòng)干預(yù)后, 甘竹、馬口潮汐動(dòng)力增強(qiáng), 潮區(qū)界、潮流界顯著上移, 這種負(fù)向波動(dòng)則逐漸消失。余水位曲率的時(shí)空演變特征如圖3c 所示?？臻g上, 曲率沿程有明顯的正負(fù)值波動(dòng), 沿程水面線呈現(xiàn)凹凸形態(tài)轉(zhuǎn)變; 而時(shí)間上, 每段曲率呈現(xiàn)波動(dòng)的階段性逆轉(zhuǎn), 與楊昊等(2020)對(duì)磨刀門河口人類活動(dòng)劃分的階段對(duì)應(yīng)較好。沿程余水位曲率的正負(fù)波動(dòng)較人類活動(dòng)前階段發(fā)生顯著變化, 三灶—燈籠山段曲率由負(fù)變正, 河道由侵蝕趨勢(shì)轉(zhuǎn)為淤積趨勢(shì); 燈籠山—竹銀段曲率由正轉(zhuǎn)負(fù), 即河道由淤積趨勢(shì)轉(zhuǎn)為侵蝕趨勢(shì); 竹銀—江門段的曲率由負(fù)轉(zhuǎn)正, 河道由侵蝕趨勢(shì)轉(zhuǎn)為淤積趨勢(shì); 江門—馬口段曲率由正轉(zhuǎn)負(fù), 河道由人類活動(dòng)前的淤積趨勢(shì)轉(zhuǎn)為侵蝕趨勢(shì)。

圖3 磨刀門河口沿程余水位Z (a)、余水位坡降S (b)和余水位曲率C (c)的演變特征圖中空白處由于1998 至2000 年江門站無水位數(shù)據(jù), 故不進(jìn)行插值。圓圈表示站位Fig. 3 Evolution of the residual water level (a), residual water level slope (b) and residual water level curvature (c) along the Modaomen estuary

為更直觀的描述強(qiáng)人類活動(dòng)前后水面坡降及曲率的季節(jié)變化及階段性差異, 本文統(tǒng)計(jì)強(qiáng)人類活動(dòng)前(1965—1990 年)與強(qiáng)人類活動(dòng)后(2001—2016 年)多年月平均水位坡降與水位曲率, 結(jié)果如圖4 所示。由圖4 可見, 夏季水位坡降大, 冬季小, 上游坡降的季節(jié)性變化大于下游, 強(qiáng)人類活動(dòng)后的調(diào)整期, 水位坡降均有不同程度的減小, 坡降的減小導(dǎo)致外海潮波傳播所受的有效摩擦減弱、河口潮差增大、咸潮上溯加劇、余流結(jié)構(gòu)改變等會(huì)進(jìn)一步導(dǎo)致曲率發(fā)生變化甚至出現(xiàn)逆轉(zhuǎn)。三灶至燈籠山年平均曲率減小0.13×10-4m·km-2, 河道傾向于全年淤積趨勢(shì); 燈籠山至竹銀年平均曲率減小1.13×10-4m·km-2, 河道傾向于夏淤冬沖趨勢(shì); 竹銀至江門段年平均曲率減小0.33×10-5m·km-2, 河道傾向于沖刷趨勢(shì); 江門至甘竹段年平均曲率增大1.30×10-4m·km-2, 河道傾向于夏沖冬淤趨勢(shì); 甘竹至馬口段年平均曲率減小1.20×10-4m·km-2, 河道傾向于夏淤冬沖趨勢(shì)。

圖4 各區(qū)段緩變期(1965—1990 年)和調(diào)整期(2001— 2016 年)內(nèi)多年月平均坡降S (a～e)與曲率C (f～j)變化Fig. 4 Monthly averaged variations in residual water level slope (a～e) and curvature (f～j) in the pre-human period(1965-1990) and post-human period (2001-2016)

4 水位空間分布影響因素

水位空間結(jié)構(gòu)呈現(xiàn)階段性變化, 主要由地形邊界影響和動(dòng)力邊界影響共同作用, 本文采用上下動(dòng)力邊界驅(qū)動(dòng)的線性回歸模型, 計(jì)算得出沿程水位的預(yù)報(bào)結(jié)果, 從而定量辨識(shí)地形和動(dòng)力邊界對(duì)水位的影響程度, 進(jìn)一步辨識(shí)出水位坡降和曲率的影響程度。模型率定計(jì)算結(jié)果(表1)和調(diào)整階段預(yù)測(cè)結(jié)果(圖 5)表明, 回歸模型率定結(jié)果良好,緩變階段(1965—1990 年)和調(diào)整階段(2001—2016年)各站位RMSE 均在0.03～0.38m,R2均大于0.89,表明該模型能較好地模擬磨刀門河口的水位時(shí)空變化, 其誤差主要取決于回歸模型建立時(shí)的誤差,如表1 所示, 誤差范圍在6%～11%, 江門站誤差最大。利用緩變階段的回歸模型系數(shù), 采用 2001—2016 年上下游動(dòng)力邊界來預(yù)估河道地形條件不變情況下的水位, 并計(jì)算實(shí)測(cè)與預(yù)測(cè)之差值, 此差值即為地形邊界變化引起的水位變化量, 從圖 5中可以看出, 2001—2016 年, 回歸模型預(yù)測(cè)水位與實(shí)測(cè)水位有較大差異, 洪水期時(shí), 其差值在0.2～2m, 遠(yuǎn)超出誤差范圍, 表明強(qiáng)人類活動(dòng)對(duì)磨刀門河口水位變化有顯著影響。

表1 磨刀門河口雙變量三參數(shù)線性回歸模型的率定結(jié)果Tab. 1 Calibration of the bivariate linear regression model with three parameters in the Modaomen estuary

圖5 磨刀門河口緩變階段回歸模型的率定(a、c、e、g、i)和調(diào)整階段水位變化過程的反演(b、d、f、h、j)a、b 為馬口站; c、d 為甘竹站; e、f 為江門站; g、h 為竹銀站; i、j 為燈籠山站。圖中虛線為1:1 率定中線Fig. 5 Calibration of the regression model in the pre-human period (a, c, e, g, i) and reconstruction of the water level in the post-human period (b, d, f, h, j) in the Modaomen estuary

使用式(6)～(8)分別計(jì)算出水位、水位坡降、水位曲率的 ΔTOT、ΔGEO、ΔBOU沿程變化, 結(jié)果如圖6和表2 所示。由圖6a 可見, 強(qiáng)人類活動(dòng)后, 全區(qū)段水位下降, 且上游較下游下降程度更明顯, 上游流量輸入減小和底床下切的疊加效應(yīng)導(dǎo)致整體水位下降, 而下游燈籠山站挖沙活動(dòng)引起的水位下降效應(yīng)與口門圍墾、海平面變化引起的水位抬升效應(yīng)相抵消。不同程度的水位下降將反映在水位坡降的沿程變化上。由圖6b 可見, 強(qiáng)人類活動(dòng)后, 全區(qū)段水位坡降減小。三灶至竹銀段水位坡降略有減小, 平均坡降減小0.22×10-2m·km-1; 竹銀站附近及河口上游人為采沙劇烈, 導(dǎo)致水位坡降急速減小, 竹銀至江門、甘竹至馬口段平均水位坡降減小 0.9 6×10-2m·km-1; 江門至甘竹段坡降略有增大, 平均坡降增大0.001×10-2m·km-1, 整體上地形邊界影響量ΔGEO與總變化值 ΔTOT趨勢(shì)一致, 動(dòng)力邊界影響量ΔBOU沿程變化不明顯, 對(duì)坡降的影響為負(fù), 僅在口門處占總變化的比值較高。圖6c 顯示, 強(qiáng)人類活動(dòng)后, 曲率的變化較為復(fù)雜, 口門段(三灶至竹銀)和河口上段(甘竹至馬口)水位曲率減小, 分別減小0.41×10-4m·km-2和1.04×10-4m·km-2; 河口中段(竹銀至甘竹) 水位曲率增大, 平均曲率增大0.21×10-4m·km-2。地形邊界影響量 ΔGEO與總變化值ΔTOT基本一致, 動(dòng)力邊界對(duì)曲率基本無影響, 表明水位曲率變化主要由強(qiáng)人類活動(dòng)導(dǎo)致的河床地形變化引起, 而圍墾、挖沙、疏浚等人類活動(dòng)導(dǎo)致的水面坡降減緩、曲率正負(fù)轉(zhuǎn)變是不可逆的。上述結(jié)果表明調(diào)整階段后珠江磨刀門河口的沖淤演變、潮波傳播等已經(jīng)處于新的格局。

圖6 磨刀門河口水位Z (a)、水位坡降S (b)、水位曲率C (c)變化值的沿程變化Fig. 6 Longitudinal variations in the residual water level (a), residual water level slope (b) and the residual water level curvature (c) along the Modaomen estuary

表2 強(qiáng)人類活動(dòng)對(duì)沿程水位、坡降、曲率的影響程度Tab. 2 Influence of strong human interventions on the residual water level, slope and curvature along the estuary

通過前文對(duì)水位空間分布的分析, 可知徑流向下游傳播過程中存在一個(gè)坡降變緩的區(qū)域, 枯季時(shí)甚至出現(xiàn)倒坡降的現(xiàn)象。撇開地形變化對(duì)水面形態(tài)的影響, 探究?jī)H由流量、外海海平面(以三灶站月平均水位作為代表)兩個(gè)因子驅(qū)動(dòng)對(duì)水位空間分布的影響。由圖7 可見, 流量對(duì)上游水位的影響較大, 而外海海平面僅對(duì)下游附近水位的影響較大。正是由于流量、海平面高度與沿程水位的非線性耦合關(guān)系,引起坡降在空間上的非線性變化, 即從口門至上游,水位坡降沿程增大, 在距口門 80～100km 處減緩,100km 處后再增大, 其中流量越大, 水位坡降也隨之增大。在下游水位項(xiàng)設(shè)置為0m 的條件下, 當(dāng)流量低于4500m3·s–1時(shí), 在中部河段(江門—甘竹)甚至?xí)霈F(xiàn)坡降為負(fù)的情況, 這是由于河流下泄受到潮流的頂托作用, 形成回水帶, 表現(xiàn)為局部壅高現(xiàn)象,故坡降在此處會(huì)有所減緩。特別是在低流量條件下,全河道整體水位均較低, 潮汐動(dòng)力相對(duì)強(qiáng)一些, 使得部分下游水位高于上游水位, 故形成倒坡降。坡降的沿程分布反映到曲率的變化上, 從口門至上游,曲率沿程變化為正-負(fù)-正, 隨著流量的增大和海平面的上升, 曲率均沒有發(fā)生逆轉(zhuǎn), 表明上下游動(dòng)力邊界驅(qū)動(dòng)無法引起曲率的正負(fù)逆轉(zhuǎn)。但在高強(qiáng)度的人類活動(dòng)后, 水面坡降減緩, 潮汐動(dòng)力增強(qiáng), 其負(fù)向波動(dòng)消失, 曲率出現(xiàn)正負(fù)逆轉(zhuǎn), 已然超出上下游動(dòng)力邊界對(duì)水面形態(tài)的驅(qū)動(dòng), 說明曲率的正負(fù)逆轉(zhuǎn)是由地形變化反饋到動(dòng)力結(jié)構(gòu)上導(dǎo)致的。

圖7 水位Z (a、b)、坡降S (c、d)、曲率C (e、f)隨流量和下游水位變化的等值線圖圖c 中紅色實(shí)線為對(duì)應(yīng)坡降為0 的等值線Fig. 7 Contour plot of residual water level (a, b), slope (c, d) and curvature (e, f) as a function of river discharge and downstream water level. The red line is the contour with zero residual water level slope

5 討論

基于韓志遠(yuǎn)等(2010)、劉鋒等(2011)、張子昊(2018)、Liu 等(2019)和Yan 等(2020)對(duì)磨刀門河口地形演變的統(tǒng)計(jì)結(jié)果, 歸納曲率的正負(fù)與河道沖淤的對(duì)應(yīng)關(guān)系, 結(jié)果如表3 所示。其中, 平均水深變化為正表示水深增大、河床呈沖刷趨勢(shì); 反之, 平均水深變化為負(fù)表示水深減小、河床呈淤積趨勢(shì)。1962至1977 年間, 河道呈現(xiàn)中上段沖刷、口門段淤積,其中燈籠山至百頃頭平均水深減小, 河槽容積因淤積減小15.2×106m3, 此河段C>0; 百頃頭至馬口段平均水深增大, 河槽容積因沖刷增大33.7×106m3,水面曲率C<0。從水位反演結(jié)果來看, 中下段通過余水位曲率可以較好的反演其沖淤趨勢(shì)。1977 至1999 年間, 河道沖淤趨勢(shì)較1977 年前發(fā)生空間上的逆轉(zhuǎn), 中下段及口門段(燈籠山至百頃頭)受人類活動(dòng)影響轉(zhuǎn)為沖刷, 但水位曲率的調(diào)整并沒有及時(shí)反饋, 此時(shí)水位曲率C>0, 指示河道仍呈現(xiàn)淤積趨勢(shì)。1999 至2005 年間, 河道呈現(xiàn)中上段沖刷、口門段淤積, 其中燈籠山至竹銀段平均水深減小0.06m,竹銀至馬口段平均水深增大0.59m; 全河道中下段河床演變趨勢(shì)與曲率對(duì)應(yīng)良好。2005 至2013 年間,河道呈現(xiàn)中段淤積、上下段沖刷, 其中百頃頭至甘竹段平均水深減小0.02m, 其他河段平均水深均增大, 全段沖淤趨勢(shì)與曲率反演結(jié)果對(duì)應(yīng)良好。綜上可見, 實(shí)測(cè)河道不同段的平均沖淤趨勢(shì)與水位曲率變化基本一致。與曲率對(duì)應(yīng)較差的河段基本是曲率為正, 而實(shí)測(cè)地形顯示為沖刷趨勢(shì), 平均沖淤趨勢(shì)與水位曲率對(duì)應(yīng)不上的原因可能是: 1) 河床的沖淤受多種因素影響, 人類高強(qiáng)度挖沙活動(dòng)使得整體河床向窄深化發(fā)展, 人為活動(dòng)對(duì)河床的不均勻改造比水位反演來的更直接, 且水位反演有一定的滯后效應(yīng); 2) 公式(13)將平均水深作為檢測(cè)沖淤變化的標(biāo)準(zhǔn), 試想, 河寬變窄平均水深變深, 而過水?dāng)嗝婷娣e變小的情況下, 實(shí)質(zhì)上是淤積的情況, 可是平均水深變深被判定為沖刷, 沖刷和淤積還是不能僅單一的從平均水深變化來進(jìn)行判斷, 需結(jié)合過水?dāng)嗝婷娣e; 3) 公式(13)的推導(dǎo)基于均勻流和平衡輸沙假設(shè)得到, 而實(shí)際輸沙并不一定是處于平衡輸沙, 比如珠江洪水期來沙量變化、人類活動(dòng)挖沙引起水體含沙量變化、地形不均勻下切引起的局部環(huán)流等因素導(dǎo)致河道非平衡輸沙, 會(huì)對(duì)沖淤判斷造成一定的誤差, 在實(shí)際應(yīng)用中有一定的局限性; 同時(shí)人類活動(dòng)河道挖沙所形成的深坑, 盡管水深變大, 但坑內(nèi)的水流是趨于靜止, 余水位曲率不會(huì)有所影響。

表3 磨刀門河口不同河段曲率與沖淤趨勢(shì)的對(duì)應(yīng)關(guān)系Tab. 3 Correlation between curvature and tendency of erosion and deposition in different reaches of the Modaomen estuary

20 世紀(jì)60 年代中期起, 磨刀門口門地區(qū)先后建設(shè)白藤堵海大堤、橫洲水道東西導(dǎo)堤和洪灣水道南北導(dǎo)堤, 口門區(qū)大片水域圍墾成陸, 形成橫洲水道和洪灣水道“一主一支”的入海格局, 加之聯(lián)圍筑閘和航道整治, 使得主流集中、流速增大、水流挾沙能力提高, 導(dǎo)致河床產(chǎn)生沖刷。20 世紀(jì)80 年代至21 世紀(jì)初, 大規(guī)模且高強(qiáng)度的河道采沙導(dǎo)致河道逐漸向窄深化發(fā)展。地形變化反饋到河道整體動(dòng)力結(jié)構(gòu)上, 口門圍墾導(dǎo)致下游過水?dāng)嗝孀冋? 且口門外海平面呈上升趨勢(shì), 水位空間分布為適應(yīng)縮窄斷面調(diào)整, 形成上凸形態(tài), 下游曲率轉(zhuǎn)為負(fù)值, 此空間動(dòng)力結(jié)構(gòu)的調(diào)整又驅(qū)使下游沖淤趨勢(shì)發(fā)生逆轉(zhuǎn), 口門段由原來的淤積轉(zhuǎn)為沖刷; 采沙使整體河床下切,航道整治擴(kuò)寬過水?dāng)嗝? 水位空間分布為適應(yīng)地形調(diào)整, 竹銀至江門段水位坡降減小0.88×10-2m·km-1,而江門至甘竹段水位坡降增大0.17×10-2m·km-1, 形成江門下凹甘竹上凸形態(tài), 高阻點(diǎn)向上游遷移, 中間段曲率轉(zhuǎn)為正值, 上游段曲率轉(zhuǎn)為負(fù)值, 此空間動(dòng)力結(jié)構(gòu)的調(diào)整又驅(qū)使沖淤趨勢(shì)發(fā)生逆轉(zhuǎn), 中間段由原來的沖刷轉(zhuǎn)為淤積趨勢(shì), 上段由原來的淤積轉(zhuǎn)為沖刷趨勢(shì)。人類活動(dòng)導(dǎo)致地形變化反饋到水位空間分布上(即曲率), 水位空間分布又能進(jìn)一步反演河道沖淤趨勢(shì), 基于長序列水位數(shù)據(jù), 可提供一種簡(jiǎn)易預(yù)測(cè)河床沖淤演變趨勢(shì)的方法, 為河口治理與水沙調(diào)控提供技術(shù)支撐。

6 結(jié)論

河口區(qū)域受徑流、潮汐、地形等多因素的耦合影響, 其水面線形態(tài)的時(shí)空變化復(fù)雜, 余水位曲率既是刻畫水面線沿程形態(tài)變化的參數(shù), 亦是河床地形沖淤變化的指示因子。本文基于珠江磨刀門河口6 個(gè)站點(diǎn)1965—2016 年的月均余水位數(shù)據(jù)及相應(yīng)時(shí)段馬口站的月均流量數(shù)據(jù), 通過數(shù)理統(tǒng)計(jì)方法和分段三次 Hermite 插值方法分析水位空間結(jié)構(gòu)變化特征, 采用雙變量線性回歸模型, 定量辨識(shí)地形和動(dòng)力邊界對(duì)水位坡降和曲率的影響程度, 并結(jié)合地形數(shù)據(jù), 初步從動(dòng)力角度探討潮平均態(tài)下曲率與底床沖淤趨勢(shì)的關(guān)聯(lián)性, 結(jié)果表明:

1) 在徑潮動(dòng)力影響下, 水位沿程分布具有明顯的季節(jié)性波動(dòng)。水位坡降與水位曲率洪枯季差異明顯, 且枯季在江門—甘竹段水位坡降存在負(fù)向波動(dòng)。強(qiáng)人類活動(dòng)影響下余水位沿程分布呈現(xiàn)階段性異變特征, 其中上段水位顯著下降, 水面線坡降變緩, 江門—甘竹段負(fù)向波動(dòng)逐漸消失, 各河段曲率出現(xiàn)正負(fù)逆轉(zhuǎn), 即從口門至上游, 水面線形態(tài)由下凹(C>0)-上凸(C<0)-下凹轉(zhuǎn)變?yōu)樯贤?下凹-上凸,河口沖淤趨勢(shì)也隨之改變。

2) 強(qiáng)人類活動(dòng)引起河床地形下切是水面空間形態(tài)發(fā)生逆轉(zhuǎn)的主要原因。河床下切引起余水位不均勻下降0.015～0.450m; 引起中上段(竹銀至馬口段)和下段( 三灶至竹銀段) 余水位坡降分別下降0.46×10-2m·km-1和0.128×10-2m·km-1; 口門段(三灶至竹銀段)和河口上段(甘竹至馬口段)水位曲率分別減小0.41×10-4m·km-2和1.04×10-4m·km-2, 河口中段(竹銀至甘竹段)水位曲率增大0.21×10-4m·km-2。上游流量和海平面高度的驅(qū)動(dòng)僅能改變整體坡降的大小, 而地形異變是導(dǎo)致坡降沿程變化趨勢(shì)改變的主要原因。

3) 水面線的凹凸性質(zhì)(即曲率的正負(fù))能較好地指示河床沖淤變化趨勢(shì), 河段平均曲率為正, 指示河床趨于淤積, 曲率為負(fù), 河床趨于沖刷侵蝕。珠江磨刀門河口大規(guī)模且高強(qiáng)度人為活動(dòng)導(dǎo)致的河床地形下切, 驅(qū)使水面坡降不均勻變緩, 水面曲率空間上發(fā)生逆轉(zhuǎn); 動(dòng)力結(jié)構(gòu)的調(diào)整驅(qū)使沖淤趨勢(shì)發(fā)生逆轉(zhuǎn), 下段和上段由原來的淤積轉(zhuǎn)為沖刷趨勢(shì), 中間段由原來的沖刷轉(zhuǎn)為淤積趨勢(shì)?；陂L序列水位數(shù)據(jù), 分析河口沿程水位分布特征及其階段性變化,可探究人類活動(dòng)、動(dòng)力結(jié)構(gòu)、地形變化之間的耦合反饋過程, 并提供一種簡(jiǎn)易預(yù)測(cè)河床沖淤演變趨勢(shì)的方法, 為河口治理與水沙調(diào)控提供技術(shù)支撐。

需要指出的是, 本文發(fā)現(xiàn)強(qiáng)人類活動(dòng)前枯季江門段月平均水位出現(xiàn)逆坡降, 但由于影響因素(比如潮非線性相互作用、江門水閘、冬季風(fēng))的不確定性,未對(duì)其形成機(jī)理進(jìn)行深入分析, 這是本文的不足之處, 也是接下來需要通過珠江河網(wǎng)數(shù)值模擬來解決的重點(diǎn)問題。