• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    珠江磨刀門河口潮波振幅梯度與上下游動力邊界的關(guān)系異變研究*

    2022-05-30 01:15:40歐素英蔡華陽楊清書
    海洋與湖沼 2022年3期
    關(guān)鍵詞:分潮馬口河口

    李 博 楊 昊 歐素英 蔡華陽①劉 鋒 楊清書

    (1. 中山大學(xué)海洋工程與技術(shù)學(xué)院 河口海岸研究所 廣東廣州 510275; 2. 河口水利技術(shù)國家地方聯(lián)合工程實驗室 廣東廣州510275; 3. 廣東省海岸與島礁工程技術(shù)研究中心 廣東廣州 510275; 4. 南方海洋科學(xué)與工程廣東省實驗室(珠海) 廣東珠海519000)

    河口是陸海徑潮動力耦合的獨特區(qū)域, 既有上游徑流匯入, 又有外海潮波上溯, 還受人類活動疊加效應(yīng)的影響, 導(dǎo)致其動力結(jié)構(gòu)復(fù)雜多變。河口動力的典型特征包括徑潮相互作用、斜壓效應(yīng)、波流耦合等,其中, 徑潮相互作用對洪水下泄、鹽淡水混合、最大渾濁帶形成、營養(yǎng)鹽和污染物輸移等物理、化學(xué)和生物過程均有直接影響。潮波振幅梯度, 即潮波振幅沿程的平均變化率, 綜合反映了徑流、潮汐和底床摩擦等多因子影響在河口空間上的差異, 是探討河口復(fù)雜徑潮動力結(jié)構(gòu)的有效切入點。因此, 揭示潮波振幅梯度與上下游動力邊界(即上游徑流和外海潮汐)的耦合過程及機制是河口動力學(xué)研究的基礎(chǔ)科學(xué)前沿問題, 可為河口海岸區(qū)域自然災(zāi)害(如咸潮上溯、洪澇、赤潮等)的防治和水資源的高效開發(fā)利用等提供重要科學(xué)依據(jù)(Schuttelaarset al, 2013; Caoet al, 2020; 楊昊等, 2020)。

    國內(nèi)外針對河口徑潮相互作用的研究由來已久。研究表明, 當潮波受到徑流和地形的調(diào)制作用時呈現(xiàn)出明顯的非線性變化(Hoitinket al, 2017; Jiet al,2019), 河口潮波衰減主要通過流量和底床的摩擦效應(yīng)實現(xiàn)(Godin, 1985; Savenije, 1992; Horrevoetset al,2004)。當流量增大時, 全日分潮和半日分潮的振幅以不同幅度衰減, 而河道底床下切, 水深增大可導(dǎo)致底床摩擦效應(yīng)減弱, 引起沿程潮波振幅不同程度的增加(Caiet al, 2012)。對河口地形進行概化并考慮下泄徑流影響, 采用一維潮波傳播解析模型可通過潮波振幅梯度和余水位梯度等參數(shù)揭示潮波傳播的時空變化特征及其動力學(xué)機制(Caiet al, 2016)。當河口動力以徑流作用為主時, 傳統(tǒng)調(diào)和分析(Pawlowiczet al,2002)的水位預(yù)報誤差明顯增大。針對流量影響下河口潮波傳播的非定常和非線性問題, 對水位序列采用連續(xù)小波變換(Guoet al, 2015; Moftakhariet al,2016)、經(jīng)驗?zāi)B(tài)分解(Panet al, 2018a)、經(jīng)驗正交函數(shù)(Panet al, 2019)、非平穩(wěn)潮汐調(diào)和分析(Matteet al,2013, 2014; Panet al, 2018b)等處理方法探討流量對潮波傳播的影響過程及機制均取得較好效果。

    近年來, 流量對潮波傳播的衰減及其閾值效應(yīng)研究得到廣泛關(guān)注。研究表明, 在長江河口感潮河段,流量與潮波振幅梯度之間存在閾值效應(yīng), 即流量小于臨界閾值時, 隨著流量增大潮波衰減效應(yīng)增強, 但當潮波振幅梯度達到極小值后, 隨著流量增大衰減效應(yīng)反而減弱(Caiet al, 2019; 張先毅等, 2019)。在珠江磨刀門河口也觀察到類似現(xiàn)象, 即在甘竹-馬口河段潮波振幅梯度隨流量變化有先減小后上升的特點,在人類活動干預(yù)較弱的自然演變階段, 潮波振幅梯度在流量為5 000~10 000 m3/s 時達到極小值, 而在人類活動強烈影響后的恢復(fù)調(diào)整階段, 該流量閾值上升至15 000~20 000 m3/s, 且其對應(yīng)的潮波振幅梯度極小值增大(Yanget al, 2020; 張先毅等, 2020), 反映因采砂導(dǎo)致的河槽容積增大對潮波振幅梯度閾值的增大效應(yīng)。綜上所述, 對于河口徑潮相互作用及潮波振幅梯度與流量動力邊界關(guān)系的研究已有較多成果,但是對于定量辨識流量驅(qū)動下潮汐調(diào)和分析結(jié)果的階段性演變以及潮波振幅梯度與下游動力邊界關(guān)系的研究等相關(guān)科學(xué)問題還尚待深入。

    位于我國粵港澳大灣區(qū)國家戰(zhàn)略核心經(jīng)濟圈的珠江三角洲是全球受到人類活動影響程度最大的區(qū)域之一。高強度的人類活動(如水庫建設(shè)、口門圍墾、無序采砂、航道疏浚等)導(dǎo)致來水來沙和地形地貌發(fā)生明顯變化, 進而顯著改變珠江河網(wǎng)的徑潮動力格局。水庫建設(shè)使河流挾沙量銳減, 三角洲河道下切加劇; 地形變化導(dǎo)致西江、北江的水沙分配格局發(fā)生變化, 河網(wǎng)中上部地區(qū)余水位梯度變緩(Liuet al, 2019);口門圍墾、口門整治疊加海平面上升的影響則導(dǎo)致磨刀門河口的形態(tài)向窄深化趨勢演變(Tanet al, 2019)。以上變化導(dǎo)致珠江主要出??? 即磨刀門河口的潮流界和潮區(qū)界顯著向上游移動, 咸潮上溯災(zāi)害風(fēng)險增大, 對粵港澳大灣區(qū)城市的居民生活、工農(nóng)業(yè)生產(chǎn)用水的水質(zhì)安全形成嚴峻挑戰(zhàn)。因此, 本文針對強人類活動驅(qū)動下的徑潮動力異變問題, 采用數(shù)據(jù)驅(qū)動模型探究流量驅(qū)動下珠江磨刀門河口潮波傳播特征的異變過程, 并分析潮波調(diào)和常數(shù)的階段性演變, 揭示潮波振幅梯度對上下游動力邊界的響應(yīng)過程及機制, 研究成果可為受到高強度人類活動影響的河口三角洲的治理和保護等提供科學(xué)依據(jù)。

    1 研究區(qū)域概況

    珠江河網(wǎng)(圖1)流域面積達4.5×105km2, 具有三江(東江、西江、北江)匯流, 八口(自東向西依次為虎門、蕉門、洪奇門、橫門、磨刀門、雞啼門、虎跳門、崖門)入海的特點, 河網(wǎng)縱橫交錯導(dǎo)致徑潮動力相互作用過程復(fù)雜。珠江每年約有2 823×108m3淡水流量及72.4×106t 懸沙量流入南海(Liuet al, 2018), 平均流量約為8 952 m3/s, 口門附近潮差為1.0~1.7 m。磨刀門河口為西江入??? 是八大口門中輸水輸沙量最大的河口, 該河口屬典型的河控型河口。據(jù)三角洲頂端馬口水文站1959~2016 年月均流量資料, 其多年平均流量為7 078 m3/s, 月均流量最大值為29 000 m3/s(6 月), 最小值為1 210 m3/s (1 月), 洪枯季差異明顯(楊昊等, 2020)。潮汐以不規(guī)則半日潮為主, 日不等現(xiàn)象顯著, 燈籠山站平均漲落潮歷時分別為5.37 h 和7.25 h, 口門處三灶站潮型系數(shù)為1.30。

    圖1 研究區(qū)域及所用潮位站或水文站位置Fig.1 The study area and the deployment of the tidal gauging stations and the hydrological station

    隨著經(jīng)濟社會高速發(fā)展, 強人類活動已成為河口三角洲系統(tǒng)演變的第三驅(qū)動力(陳吉余等, 2008)。口門圍墾、聯(lián)圍筑閘、河道采砂及上游流域水庫建設(shè)等高強度人類活動極大影響磨刀門河口的河道形態(tài)和徑潮動力, 其“動力-沉積-地貌”體系發(fā)生顯著變化(吳超羽, 2018)。20 世紀80 年代前, 磨刀門地區(qū)主要人類活動為大規(guī)模灘涂圍墾, 至 1 9 9 7 年已有16.5 km2的淺灘被圍填成陸地(Tanet al, 2019)。圍墾、疏浚、口門整治等一系列工程措施使河口口門區(qū)水域面積減小, 入海水道水深增大, 至20 世紀90 年代初入??陂T向海延伸約13 km (何用等, 2018)。其次, 20世紀50~70 年代中期, 大規(guī)模的聯(lián)圍筑閘導(dǎo)致三角洲腹部水位壅高, 中上游水面比降變緩, 流速減小, 水流挾沙能力降低, 河床產(chǎn)生淤積(劉鋒等, 2011)。20世紀80 年代中期開始, 由于用沙需求激增, 在珠江河網(wǎng)區(qū)出現(xiàn)大規(guī)模河道采砂活動。磨刀門采砂活動主要集中在珠海大橋以北河段, 1991~2000 年共挖沙420×104m3(韓志遠等, 2010)。截至2006 年, 西江流域大壩總庫容達5.38×1010m3, 占其入海淡水流量的24.7% (Liuet al, 2018)。大壩建成后流域內(nèi)大部分沉積物被截留在上游地區(qū), 西江干支流輸沙量均呈下降趨勢, 下游泥沙供給量銳減(Wanget al, 2011)。在上述強人類活動的綜合影響下, 珠江磨刀門河口沿程河床大幅下切, 河槽容積增大, 導(dǎo)致近年來潮汐動力顯著增強, 咸潮上溯加劇(韓志遠等, 2010; 張先毅等, 2020)。

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

    2.1 數(shù)據(jù)及其處理方法

    本文所用數(shù)據(jù)取自《廣東省水文年鑒》第八卷和廣東省水文局, 潮位原始數(shù)據(jù)高程為凍結(jié)基面, 已統(tǒng)一轉(zhuǎn)換至珠江基面。數(shù)據(jù)包括磨刀門河口沿程4 個潮位站(三灶、燈籠山、竹銀、甘竹站)逐日高、低潮位數(shù)據(jù)和相應(yīng)時段馬口水文站日均流量及逐日高、低潮位數(shù)據(jù), 并對潮位和流量數(shù)據(jù)進行分段三次Hermite插值處理后獲得的逐時數(shù)據(jù)。站點具體信息見表1。

    表1 磨刀門河口各站點數(shù)據(jù)信息Tab.1 Information of the adopted tidal gauging stations and hydrological station in the Modaomen estuary

    2.2 R_TIDE 數(shù)據(jù)驅(qū)動模型

    本文采用徑潮耦合的數(shù)據(jù)驅(qū)動模型(river-tidal harmonic analysis, R_TIDE)對逐時潮位數(shù)據(jù)進行調(diào)和分析(歐素英等, 2017), 模型可輸出重構(gòu)的逐時水位、特定分潮的振幅和遲角等。在徑流動力占主導(dǎo)的河口,潮波上溯受到徑流季節(jié)性周期變化的調(diào)制影響, 傳統(tǒng)的潮汐調(diào)和分析模型對河口區(qū)潮波特征的分析和預(yù)報誤差較大, 不能滿足徑潮信號分離的要求。為有效辨識潮波在流量影響下的傳播特征, 基于Matte 等(2013, 2014)提出的非平穩(wěn)調(diào)和分析(nonstationary tidal harmonic analysis, NS_TIDE)思路, 假設(shè)河口任意位置潮波振幅和遲角主要受上游流量和地形變化的非線性調(diào)制影響, 在此基礎(chǔ)上進一步假定分析時段內(nèi)河口地形邊界不變, 將信號分成由流量Q引起的水位變化和徑流調(diào)制影響下的潮水位變化, 通過信號分離定量分析感潮河段徑潮動力相互作用。

    無徑流調(diào)制時, 對于月球和太陽引起的周期性潮汐現(xiàn)象可看作是許多假想天體引起的簡單的潮汐簡諧波動的總和(Pawlowiczet al, 2002), 即潮位曲線可近似看作是由許多振幅、周期、相位不同的分潮組成。潮位可表達成以下多個余弦函數(shù)的疊加形式:

    其中,z(t)為站點的實測潮位;t為時間;z0為平均海面高度;N為分潮個數(shù);fi為分潮振幅的訂正因子, 為時間的函數(shù), 常取一年的中值;Hi為分潮平均振幅;σi為分潮角頻率;Vi為分潮初相角;ui為天文相角的交角訂正角,gi為由于海底摩擦、海水慣性引起的遲角。其中Hi和gi合稱為分潮調(diào)和常數(shù)。

    Kukulka 等(2003)基于一維圣維南方程組推導(dǎo)得出河口三角洲內(nèi)任意位置x的潮波振幅變化:

    其中,j=1,2,3,…,m,m為逐時序列的數(shù)據(jù)個數(shù);r為衰減系數(shù);ε0(tj)為河口口門位置x=0 處的潮波振幅;a'為常數(shù)。衰減系數(shù)r為潮波傳播速度、水深、流量和河道摩擦系數(shù)的函數(shù)(Jayet al, 1997)。若摩擦系數(shù)采用曼寧公式, 且摩擦系數(shù)與潮波傳播速度均為水深的函數(shù), 則衰減系數(shù)r可簡寫為水深h(tj)和流量Q(tj)的函數(shù):

    或:

    其中,h(tj)和Q(tj)表示隨時間變化的水深和流量;p0、p1、p2、β、γ均為待定參數(shù)。式(3)未考慮外海非潮動力邊界比如風(fēng)暴等因素的影響, 這是與Matte 等(2013,2014)提出的非定常調(diào)和分析模型的不同之處。假定研究時間段河道底床高程不隨時間變化, 則式(4)中水深項可歸到常數(shù)項, 衰減系數(shù)r可簡寫為

    即河口潮汐非線性變化僅受流量驅(qū)動影響。結(jié)合式(1)、(2)、(5)可得流量影響下任意位置x的潮水位變化:

    式(6)的待定系數(shù)bi、ci、γ可通過粒子群優(yōu)化算法(Kennedyet al, 1995)或求解帶約束的非線性多變量方程組來確定該參數(shù)的最優(yōu)值。由于流量對不同分潮的影響程度存在差異(Guoet al, 2015), 半日分潮的衰減效應(yīng)大于全日分潮, 將流量對每一個分潮族的調(diào)制影響均采用不同的待定系數(shù), 則式(6)可表示為:

    其中,

    式(7)即為流量驅(qū)動下的R_TIDE 模型,S(tj)為底床高程或水深變化、海平面、流量等引起的平均水面變化, 稱為余水位項;F(tj)代表不同流量條件下徑潮耦合引起的水位變化, 稱為徑潮耦合項。d0、d1、v0、v1、r0、r1為模型的待定參數(shù), 通過求解帶約束的非線性多變量方程組來確定。同時, 由于磨刀門河口徑流動力占優(yōu)勢導(dǎo)致徑潮動力非線性因素作用較強, 洪水流量時導(dǎo)致三角洲上段為徑流控制, 無潮汐波動。因此, 在R_TIDE 模型中引入臨界流量QC, 當連續(xù)2 d 潮差小于某個值(如0.001 m)時所對應(yīng)的最小流量即為臨界流量, 當流量大于QC時, 表示該站點及其以上河段不存在潮汐信號。

    假設(shè)河口三角洲任意位置采用時間間隔為Δt的觀測潮位y(tj)(j=1,2,3,…,m)序列, 根據(jù)實測資料序列長度mΔt、分潮頻率差和 Rayleigh 準則判據(jù) Δσ=max[(mΔt)-1,σi-σj], 選擇待提取的分潮, 將參數(shù)tj、σj、Q(tj)代入式(9), 要求模型計算結(jié)果y(tj)和實測水位z(tj)的誤差平方和達到最小, 即:

    選擇對應(yīng)資料長度的全日分潮、半日分潮等分潮族, 對實測潮位進行回歸擬合, 回歸模型效果以均方根誤差ERMS和相關(guān)指數(shù)R2表示, 即:

    2.3 徑潮動力的特征指標

    3 結(jié)果

    3.1 徑潮動力格局的階段性演變

    前期研究表明, 基于燈籠山-馬口河段月均余水位梯度(或潮波振幅梯度絕對值)與流量的雙累積曲線變化趨勢, 可知曲線斜率在1990 年與2000 年發(fā)生明顯改變(張先毅等, 2020)。根據(jù)上述結(jié)果, 可將磨刀門河口徑潮動力格局的演變分為三個階段: 1990 年前為自然演變和人類活動共同作用的階段, 但此階段人類活動尚未對河口徑潮動力格局產(chǎn)生明顯影響, 河口處在動態(tài)平衡期, 以自然演變?yōu)橹? 稱為“平衡期”(簡稱P1); 2000 年后, 河道中下游采砂活動基本已停止, 灘涂圍墾面積銳減, 河口在人類活動影響下已發(fā)生徑潮動力格局轉(zhuǎn)換, 為強人類活動干預(yù)后河口的自適應(yīng)調(diào)整階段, 稱為“調(diào)整期”(簡稱P2); 1991~2000年為人類活動影響最強烈的時期, 河口從自然演變?yōu)橹鞯钠胶鈶B(tài)向強人類活動干預(yù)后的平衡態(tài)轉(zhuǎn)變,為兩個平衡態(tài)的過渡階段, 稱為“過渡期”。

    3.2 R_TIDE 數(shù)據(jù)驅(qū)動模型率定效果

    圍繞珠江磨刀門河口徑潮動力異變問題, 采用流量驅(qū)動的R_TIDE 數(shù)據(jù)驅(qū)動模型對徑潮信號進行分離。表2 為磨刀門河口兩個階段沿程站點的臨界流量QC、模型主要率定參數(shù)及相關(guān)的模型評價指標(ERMS和R2)結(jié)果。C1~C9為不同分潮族所使用的參數(shù), 其中C1代表長周期氣象分潮, 其周期約為15 d, 反映一個大小潮周期變化;C2、C3分別代表全日和半日分潮族,為主導(dǎo)該區(qū)域潮汐動力的分潮族;C4~C9分別代表短周期波動, 由分潮淺水變形效應(yīng)(如M4、MS4、M6、M8)產(chǎn)生, 反映地形、底床摩擦等邊界條件的影響。由表2 可知, 越往下游ERMS值越大, 表明模型對潮流優(yōu)勢區(qū)域水位的重構(gòu)效果略低于上游區(qū)域, 但是其R2均大于0.79, 模型結(jié)果合理。對于馬口站, 由于其為流量驅(qū)動邊界, 除受洪峰影響外還受邊界效應(yīng)影響, 因此其ERMS相對較大。調(diào)整期階段的ERMS比平衡期普遍增大,R2減小, 表明強人類活動對地形的改造導(dǎo)致水面線為適應(yīng)地形變化, 水位對流量的響應(yīng)敏感度有所降低, 由此一定程度上影響了模型的準確性。

    圖2 為R_TIDE 模型在平衡期和調(diào)整期重構(gòu)得到的不同站點日均水位與實測值的對比結(jié)果, 其中黑色虛線表示實測值與重構(gòu)值完全一致。當流量大于臨界流量QC時潮波基本消失, 率定過程中流量大于QC時間段所對應(yīng)的潮波信號為虛假信號, 結(jié)果分析部分已剔除這部分信息。由圖2 可知, 馬口站下游站點模型擬合的效果較好(ERMS介于0.14~0.22 m,R2介于0.79~0.99)。但由于夏季洪峰出現(xiàn)次數(shù)多且持續(xù)時間長, 對徑潮信號分離產(chǎn)生較大影響, 模型誤差較大,因此當流量接近臨界流量QC時, 馬口站重構(gòu)水位比實測值偏大。

    圖2 R_TIDE 數(shù)據(jù)驅(qū)動模型率定結(jié)果Fig.2 Calibration of the adopted R_TIDE data-driven model

    由表2 可見, 調(diào)整期的QC和γ均大于平衡期, 主要原因為調(diào)整期階段潮汐動力增強, 潮區(qū)界向上游移動。流量指數(shù)γ反映地形邊界條件變化, 其值增大表明河床下切, 過水斷面增大。其中三灶站變化最大,反映入海水道水深劇增, 而上游各站點水深均有不同程度增大。竹銀、燈籠山站主要受采砂影響, 而馬口、甘竹站點除受采砂影響外, 還受上游流域水庫調(diào)蓄等影響, 導(dǎo)致馬口水文站的分流、分沙比發(fā)生變化, 進而顯著改變底床形態(tài)。竹銀及其下游站點為潮流優(yōu)勢河段,QC均一致, 而甘竹、馬口則受強烈徑流影響, 流量越大, 潮區(qū)界越往下游移動, 因此,QC逐漸減小。

    表2 R_TIDE 數(shù)據(jù)驅(qū)動模型的部分參數(shù)及效果評價Tab.2 Calibrated parameters derived from R_TIDE data-driven model and its model performance

    3.3 M2 分潮振幅及其梯度的階段性變化特征

    由于磨刀門水位及其梯度發(fā)生階段性異變, 導(dǎo)致潮波振幅及其梯度亦發(fā)生顯著變化。磨刀門河口潮汐動力主要由M2分潮主導(dǎo), 因此選用其作為代表性分潮進行分析。M2分潮振幅及其梯度的季節(jié)及年均變化如表3 所示。潮波向上游傳播過程中振幅沿程衰減, 洪季振幅小于枯季(三灶站除外), 表明徑流對潮波傳播存在明顯的衰減作用, 尤其是在磨刀門河口上游。平衡期燈籠山站振幅的洪枯季差異僅為0.03 m,甘竹站則為0.07 m, 調(diào)整期潮波振幅有明顯增大(三灶站除外), 且上游馬口站振幅的洪枯季差異增大(約為0.02 m), 其下游的甘竹、竹銀站則減小, 燈籠山、三灶站振幅的洪枯季差異基本不變。

    表3 磨刀門河口各潮位站M2 分潮振幅及其梯度的季節(jié)變化及年均變化Tab.3 Seasonal and annual alterations in amplitude and its gradient of M2 tidal constituent observed at each tidal gauging station along the Modaomen estuary

    為探究潮波傳播過程的時空不均勻性, 對比分析強人類活動前后各站點日均M2分潮振幅及相鄰站點間日均M2分潮潮波振幅梯度的階段性演變, 繪制時間過程線如圖3 所示, 圖中黑色實線表示階段性演變?yōu)?。圖3a 為強人類活動前后各站點日均M2分潮振幅變化值的年內(nèi)分布, 其中日均M2分潮振幅變化值定義為調(diào)整期的M2分潮振幅的日平均值減去平衡期所對應(yīng)的值, 記為Δη(以下日均M2分潮潮波振幅梯度變化值的定義類似, 記為Δδ)。由圖3a 可見, 站點距離口門越遠其值越大, 其中三灶站振幅年均減小0.01 m, 而馬口站則變化0.11 m。季節(jié)變化也有類似規(guī)律(馬口站夏季和甘竹站冬季除外), 表明振幅變化量的季節(jié)性差異有向上游累積的效應(yīng)。日均潮波振幅梯度的時空變化如圖3b 所示。潮波振幅梯度最小值出現(xiàn)在夏季, 最大值出現(xiàn)在冬季。三灶-甘竹河段的潮波振幅梯度在調(diào)整期階段的洪枯季差異減小(各河段洪枯季差異的階段性變化分別為 7.14×10-7、4.86×10-6、5.11×10-6m-1), 且越往上游減小量越大, 表明強人類活動對潮波振幅梯度的增大效應(yīng)在甘竹以下河段都有累積效應(yīng), 但在甘竹-馬口河段, 洪枯季差異雖減小但幅度小于竹銀-甘竹河段(洪枯季差異的階段性變化為4.24×10-6m-1), 表明強人類活動干預(yù)后甘竹-馬口河段潮波振幅梯度對流量響應(yīng)減弱。從潮波振幅梯度變化量來看, 三灶-燈籠山河段在秋季變化最大,燈籠山-甘竹河段則為夏季變化最大, 甘竹-馬口河段則是夏季。上述結(jié)果表明三灶-燈籠山河段水動力主要受海平面變化控制(秋季三灶余水位最大); 而甘竹-馬口河段受徑流控制為主, 由于該區(qū)域在調(diào)整期明顯受水庫調(diào)蓄影響(即洪季蓄水、枯季放水), 雖然潮波振幅梯度夏季階段性變化最大, 冬季最小, 但是其洪枯季差異減小, 是流量季節(jié)性差異減小的直接反映。

    圖3 強人類活動前后磨刀門河口日均M2 分潮潮波振幅及其梯度的階段性差異Fig.3 Stepwise alteration in tidal amplitude and gradient under intensive human activities in the Modaomen estuary

    3.4 潮波振幅梯度與上下游動力邊界的關(guān)系

    圖4 為三灶-馬口河段M2分潮日均潮波振幅梯度與馬口站日均流量的關(guān)系曲線(圖4a)及其季節(jié)變化(圖4b~4e)。平衡期階段三灶-馬口河段潮波振幅梯度極小值(記為δM)為-1.41×10-5m-1, 對應(yīng)流量閾值為10 500 m3/s, 而在調(diào)整期階段則分別變?yōu)?1.38×10-5m-1和19 500 m3/s (其中流量閾值定義為當潮波振幅梯度達到極小值所對應(yīng)的流量值, 記為QM)。日均尺度下, 閾值效應(yīng)在調(diào)整期的冬季基本消失。從季節(jié)變化來看, 平衡期的秋季潮波振幅梯度極小值最小, 調(diào)整期則為夏季最小, 冬季最大, 大流量下潮波衰減效應(yīng)更明顯。調(diào)整期階段潮波振幅梯度極小值(約為-1.38×10-5m-1)和對應(yīng)流量閾值(約為19 500 m3/s)均大于平衡期。

    圖4 三灶-馬口河段M2 分潮日均潮波振幅梯度(δ)與馬口站日均流量(Q)的閾值關(guān)系(a)及其季節(jié)變化(b~e)Fig.4 Relationship between the daily averaged M2 tidal amplitude gradient over SZ-MK reach (δ) and the river discharge (Q) observed at the MK hydrological station to illustrate the threshold effect (a) and its seasonal alteration (b~e)

    潮波振幅梯度是徑潮動力相互作用的直接結(jié)果,因此類似流量閾值效應(yīng), 外海潮汐動力變化(即周期約為15 d 的大小潮)亦會導(dǎo)致潮波振幅梯度產(chǎn)生振幅閾值效應(yīng)(記為ηM)。圖5 為磨刀門河口三灶-馬口河段日均M2分潮潮波振幅梯度與三灶站日均M2分潮振幅的關(guān)系曲線。由圖5 可見, 日均潮波振幅梯度極小值在調(diào)整期略有增大, 約為2.49×10-7m-1, 對應(yīng)的振幅閾值基本不變(約增大0.14 cm)。圖5b~5e 為日均M2分潮潮波振幅梯度與三灶站日均M2分潮振幅閾值關(guān)系的季節(jié)變化。與圖4 類似, 除調(diào)整期的冬季振幅閾值效應(yīng)消失外(圖5e), 其余季節(jié)均存在較明顯的振幅閾值效應(yīng), 其中夏季振幅梯度極小值最小, 其次為春、秋兩季, 最大為冬季。而平衡期階段秋季的振幅閾值最小, 冬季最大。調(diào)整期冬季閾值效應(yīng)消失的主要原因是下泄流量較小, 即便經(jīng)過流域水庫調(diào)蓄,調(diào)整期口門振幅依舊未能達到使衰減效應(yīng)達到最大的臨界值(約0.43 m)。達到振幅閾值前, 由于大潮流速振幅較大導(dǎo)致其底床摩擦較大, 衰減效應(yīng)增強, 潮波振幅梯度減小。而達到振幅閾值后, 流速振幅增大導(dǎo)致的潮能耗散不足以平衡地形輻聚增強(大潮時水深增大)導(dǎo)致的潮能增加, 因此, 潮波振幅梯度略微增大。圖6 為相鄰站點間沿程潮波振幅梯度隨口門三灶站振幅的變化曲線。結(jié)果表明, 振幅閾值效應(yīng)僅存在于徑流動力占優(yōu)勢的河段, 主要存在于甘竹-馬口河段(僅調(diào)整期的冬季消失), 另外, 在竹銀-甘竹河段(春、夏季和平衡期的秋季)及燈籠山-竹銀河段(調(diào)整期的夏季)也有類似現(xiàn)象, 而在三灶-燈籠山河段, 該現(xiàn)象則不存在。對于甘竹-馬口河段, 平衡期階段, 其秋季潮波振幅梯度極小值最小, 為-3.55×10-5m-1,冬季值最大, 為-3.52×10-5m-1。雖然冬季余水位減小、底床摩擦增大, 但由于流量減小, 對潮能的耗散作用亦減弱, 潮波振幅梯度反而增大, 表明在年內(nèi)變化中甘竹-馬口河段主要受徑流影響。而在調(diào)整期階段, 冬季振幅閾值效應(yīng)消失, 這是由于冬季衰減效應(yīng)減小, 尚未達到閾值(約為-3.23×10-5m-1); 其余季節(jié)均存在閾值效應(yīng), 表明要在更大的外海潮汐動力驅(qū)動下才能使該河段發(fā)生閾值效應(yīng), 這與強人類活動驅(qū)動下各站點(三灶站除外)振幅均有不同程度的抬升有關(guān), 即兩個站點間振幅差值需要達到臨界值, 潮波振幅梯度極小值才會出現(xiàn)。

    圖5 三灶-馬口河段M2 分潮日均潮波振幅梯度(δ)與三灶站日均M2 分潮振幅(ηSZ)的閾值關(guān)系(a)及其季節(jié)變化(b~e)Fig.5 Relationship between the daily averaged M2 tidal amplitude gradient (δ) in SZ-MK sector and the M2 tidal amplitude (ηSZ)observed at the SZ tidal gauging station to illustrate the threshold effect (a) and its seasonal alteration (b~e)

    圖6 相鄰站點間河段M2 分潮日均潮波振幅梯度(δ)與三灶站日均振幅(ηSZ)關(guān)系曲線的季節(jié)變化Fig.6 Seasonal alteration in the relationship between the daily averaged M2 tidal amplitude gradient of the sub-sector along neighboring tidal gauging stations and tidal amplitude observed in SZ tidal gauging station

    對于下游河段(竹銀以下), 該區(qū)域主要受口門圍墾、河道整治工程等地形邊界影響, 對徑流非線性作用的響應(yīng)不明顯, 潮波振幅梯度隨口門振幅的變化基本呈線性關(guān)系, 而燈籠山-竹銀河段出現(xiàn)閾值效應(yīng)(平衡期的夏季)主要是由于夏季衰減效應(yīng)已達到調(diào)整期階段的閾值(約為-3.92×10-5m-1), 而平衡期振幅梯度極值更小, 潮波衰減尚未達到該臨界值。竹銀-甘竹河段閾值效應(yīng)出現(xiàn)的原因與甘竹-馬口河段類似。

    表 4 為三灶-馬口河段及相鄰站點間河段的潮波振幅梯度閾值及其對應(yīng)的振幅閾值。沿程潮波振幅梯度極小值變化規(guī)律與該河段M2分潮潮波振幅梯度的空間變化規(guī)律相似, 竹銀-甘竹河段衰減效應(yīng)減弱,到甘竹-馬口河段衰減效應(yīng)又增大, 表明甘竹站以上河段受強烈徑流調(diào)制, 潮波衰減效應(yīng)增強。在上游河段振幅沿程減小, 若要達到潮波振幅梯度閾值已不需要更大的口門振幅, 因此, 對應(yīng)振幅閾值沿程略有減小。

    表4 三灶-馬口河段M2 分潮潮波振幅梯度極小值(δM)及其對應(yīng)的振幅閾值(ηM)的季節(jié)變化Tab.4 Seasonal alterations in the minimum of daily averaged M2 tidal amplitude gradient (δM) in SZ-MK sector and the critical tidal amplitude (ηM) observed at SZ tidal gauging station

    4 討論

    上述結(jié)果表明, 磨刀門河口流量-余水位及其梯度關(guān)系在調(diào)整期階段發(fā)生顯著變化, 同流量下沿程各站點余水位及其梯度明顯下降(張蔚等, 2008; Yanget al, 2020)。然而三角洲頂端的馬口水文站年均流量僅減小816 m3/s (約為11%), 下游口門處三灶站M2分潮振幅年均值僅變化0.01 m, 表明上下游動力邊界的階段性變化不是導(dǎo)致潮波振幅梯度與上下游關(guān)系異變的主要原因, 地形異變及其導(dǎo)致的底床摩擦效應(yīng)變化才是導(dǎo)致徑潮動力格局轉(zhuǎn)換的主導(dǎo)因素。

    余水位梯度作為表征徑潮動力的重要參數(shù), 其形成演變是探究地形異變對潮波傳播影響機制的重要依據(jù)和有效切入點, 其解析表達式可由一維動量守恒方程得到(Savenijeet al, 2005)

    式中,x為距離口門處的距離;ρ為水體密度;U為斷面平均流速;Z為自由水面高程;D為水深;g為重力加速度;K為曼寧摩擦系數(shù)的倒數(shù)。忽略密度效應(yīng)和對流加速度項并假定流速具有周期性, 式(15)在一個潮周期內(nèi)積分可得式(16)(Caiet al, 2014)。

    式(16)為余水位梯度的解析表達式, 其中Z為一個潮周期內(nèi)的余水位。表明余水位梯度主要與有效摩擦項相平衡(Buschmanet al, 2009; Sassiet al, 2013), 即余水位梯度越大, 潮波所受有效摩擦越大。圖7 為磨刀門河口全河段及各河段M2分潮日均潮波振幅梯度和日均余水位梯度的關(guān)系曲線。三灶-馬口河段的曲線均為下凹形, 表明余水位梯度增大過程中, 潮波振幅梯度減小速度減緩并趨近于其極小值, 達到極小值后變大, 表明潮波振幅梯度的變化已不受底床摩擦主導(dǎo)。竹銀以下河段曲線為上凸形, 表明余水位梯度增大時有效摩擦增大, 在該區(qū)域中, 靠近口門的三灶-燈籠山河段(圖7b)在調(diào)整期變化速率更快, 該變化主要由口門圍墾導(dǎo)致的水深變淺主導(dǎo), 而燈籠山-竹銀河段(圖7c)則受采砂和圍墾的共同影響。竹銀-甘竹河段(圖7d)曲線由平衡期的下凹轉(zhuǎn)變?yōu)檎{(diào)整期的上凸形, 根據(jù)該河段余水位梯度和潮波振幅梯度曲線,相同的余水位梯度條件下, 潮波振幅梯度增大, 表明有效摩擦對其影響減弱, 這與航道疏浚后水深增大有關(guān)。甘竹-馬口河段(圖7e)曲線為下凹形, 這表明影響該河段潮波傳播的多種因素耦合作用增強, 調(diào)整期余水位梯度顯著減小而水深增大, 當兩種效應(yīng)的耦合作用達到平衡時即達到潮波振幅梯度閾值, 因此調(diào)整期階段潮波振幅梯度閾值大于平衡期。

    圖7 相鄰站點間河段M2 分潮日均潮波振幅梯度(δ)與日均余水位梯度(S)的變化Fig.7 Stepwise alteration in the relationship between the daily averaged M2 tidal amplitude gradient (δ) in SZ-MK sector and the corresponding residual water level gradient (S) observed in each sub-sector

    余水位梯度的階段性變化是水面線適應(yīng)地形變化的直接反映。在人類活動干預(yù)前, 河床緩慢抬升,底坡降變化不大。1962~1977 年為緩慢淤積狀態(tài), 磨刀門河口竹洲頭-大排沙河段、大排沙-燈籠山河段平均水深分別減小0.49, 0.28 m; 1977~1999 年, 水道受到明顯沖刷, 河道深泓下切深度為0.59~1.70 m 不等,平均水深增加1.13 m (劉鋒等, 2011); 磨刀門河口沿程平均寬深比從1960 年的6.25 減小到1999 年的4.73(河道趨向窄深化趨勢發(fā)展), 同期河槽容積增大34.15×106m3(Liuet al, 2019)。表5 統(tǒng)計了三灶-竹銀河段1964、1977、1999、2016 年各水深區(qū)間的容積和面積。按水深分成0~2、2~5、5~10 和>10 m 四個區(qū)間, 分別計算各區(qū)間的河槽容積和面積。1964~1977 年, 各水深段容積和面積均顯著減小, 反映以圍墾為主導(dǎo)的淤積作用。1977~1999 年, 口門圍墾和采砂共存, 前者導(dǎo)致口門大部分水域被圍填成陸地, 0~2 m 水深區(qū)間容積和面積均銳減, 然而水深5 m 以下容積和面積顯著增大, 表明河床大幅下切,兩岸侵蝕效應(yīng)增強, 河道被沖刷, 容積和面積增加量分別為每年3.78×106m3和0.90×106m2。1999~2016年, 強人類活動逐漸減弱, 但潮波傳播過程已發(fā)生異變, 加上上游水庫建設(shè)導(dǎo)致泥沙減少, 導(dǎo)致河道水沙分配格局發(fā)生變化, 潮汐動力顯著增強, 沖刷效應(yīng)進一步加劇, 5 m 以下水深容積持續(xù)增大, 增加量為每年1.06×106m3, 面積僅增大2.93×106m2。當沖淤態(tài)勢由淤積轉(zhuǎn)變?yōu)闆_刷時, 河道容積增大, 底床摩擦減小, 潮波振幅梯度增大。

    表5 三灶-竹銀河段1964、1977、1999、2016 年各水深區(qū)間的河槽容積和面積對比Tab.5 Comparison in water volume and surface area of the channel in different water depths in SZ-ZY sector observed in 1964, 1977,1999, 2016

    5 結(jié)論

    本文基于珠江磨刀門河口沿程4 個潮位站的高、低潮數(shù)據(jù)及相應(yīng)時段馬口水文站高、低潮數(shù)據(jù)和流量數(shù)據(jù), 采用基于流量驅(qū)動的R_TIDE 數(shù)據(jù)驅(qū)動模型對河口潮波傳播的異變過程進行分析, 根據(jù)余水位及其梯度、M2分潮振幅及其梯度的階段性變化重點探討潮波振幅梯度與上下游動力邊界的關(guān)系異變過程及機制, 主要結(jié)論如下:

    (1) 與平衡期相比, 在強人類活動影響后的調(diào)整期階段, 各站點日均M2分潮振幅及其梯度增大, 其中振幅在秋季變化最大(燈籠山-馬口分別為0.06、0.10、0.10、0.13 m), 冬季最小(馬口除外, 燈籠山、竹銀、甘竹分別為0.05、0.08、0.07 m), 年內(nèi)差異不大; 潮波振幅梯度變化則較為復(fù)雜, 三灶-甘竹河段沿程減小, 甘竹-馬口河段階段性變化增大。

    (2) 日均M2分潮潮波振幅梯度與下游動力邊界之間存在閾值效應(yīng), 主要出現(xiàn)在甘竹-馬口河段。強人類活動驅(qū)動下潮波振幅梯度閾值增大(三灶-馬口河段增量約為2.50×10-7m-1), 振幅閾值基本不變(僅增大約0.14 cm)。由于底床有效摩擦減小導(dǎo)致沿程潮能耗散減弱, 因此當口門振幅相同時, 調(diào)整期的潮波振幅梯度大于平衡期, 隨著振幅增大, 底床摩擦效應(yīng)越明顯, 導(dǎo)致調(diào)整期潮波衰減速度減緩, 當兩個階段的口門振幅均達到振幅閾值時(此處假定振幅閾值相同),則調(diào)整期對應(yīng)的潮波振幅梯度大于平衡期。

    (3) 對于潮波振幅梯度與大小潮的關(guān)系異變分析需綜合考慮地形、摩擦、流量等多因素影響。由于上下游動力邊界的階段性變化不足以導(dǎo)致上述關(guān)系的異變, 因此地形及其引起的底床摩擦變化則是徑潮動力異變的根本原因。從1964~2016 年, 10 m 以下水深的容積增大15.26×106m3, 0~2 m 水深容積則減小366.94×106m3, 反映河道窄深化的演變趨勢。由于地形的異變, 在口門振幅增大過程中, 達到振幅閾值前,由于底床摩擦增大, 小潮時的潮波振幅梯度大于大潮, 而達到振幅閾值后, 此時以地形輻聚為主導(dǎo)影響因子導(dǎo)致潮波振幅梯度增大。

    猜你喜歡
    分潮馬口河口
    馬口魚
    垂釣(2023年11期)2024-01-21 16:07:04
    大亞灣雙峰水位的形成條件及準調(diào)和分量應(yīng)用的分析
    山東鄰海長周期分潮對深度基準面的影響分析
    非遺視野下湖北馬口窯的保護、繼承與開發(fā)
    流行色(2019年5期)2019-12-13 18:38:50
    來自馬口窯的對話
    ——馬口窯文獻與當代陶藝創(chuàng)作研究展
    馬口煤礦分層開采工作面上覆采空區(qū)自燃治理技術(shù)研究
    中國煤炭(2016年9期)2016-06-15 20:29:53
    他們?yōu)槭裁催x擇河口
    河口,我們的家
    特殊的河口水
    河口
    亚洲av福利一区| 91aial.com中文字幕在线观看| 欧美激情 高清一区二区三区| 午夜激情av网站| 国产一区亚洲一区在线观看| 这个男人来自地球电影免费观看 | 在线观看三级黄色| 亚洲精华国产精华液的使用体验| 国产精品久久久久久精品电影小说| 日本猛色少妇xxxxx猛交久久| 人成视频在线观看免费观看| 日本欧美视频一区| 国产乱来视频区| 亚洲色图综合在线观看| 婷婷色综合大香蕉| 久久久亚洲精品成人影院| 欧美日韩综合久久久久久| 欧美人与善性xxx| 国产精品不卡视频一区二区| 久热久热在线精品观看| 亚洲精品日本国产第一区| 国产精品三级大全| 国产极品粉嫩免费观看在线 | 91久久精品电影网| 久久鲁丝午夜福利片| 国产色爽女视频免费观看| 日本黄色片子视频| 午夜福利影视在线免费观看| 国产高清国产精品国产三级| 精品一区二区免费观看| 亚洲精品日韩在线中文字幕| 熟女电影av网| 午夜福利在线观看免费完整高清在| 熟女电影av网| 在线观看人妻少妇| 三级国产精品片| 亚洲av免费高清在线观看| 精品少妇久久久久久888优播| 三级国产精品片| 成人无遮挡网站| 97在线视频观看| 一本色道久久久久久精品综合| 欧美亚洲 丝袜 人妻 在线| 午夜福利在线观看免费完整高清在| 国产午夜精品一二区理论片| 一个人免费看片子| 嘟嘟电影网在线观看| 久久国内精品自在自线图片| 国产亚洲一区二区精品| 亚洲第一区二区三区不卡| 国产伦理片在线播放av一区| 中文字幕免费在线视频6| 在线观看美女被高潮喷水网站| 国产精品无大码| 色网站视频免费| 在线观看人妻少妇| 一边摸一边做爽爽视频免费| 日本av手机在线免费观看| 色哟哟·www| 久久久久久伊人网av| 26uuu在线亚洲综合色| 国产精品一国产av| 亚洲精品aⅴ在线观看| 国产av码专区亚洲av| 免费av不卡在线播放| 国产成人一区二区在线| 成人免费观看视频高清| 国产免费福利视频在线观看| 人成视频在线观看免费观看| 人成视频在线观看免费观看| av网站免费在线观看视频| 美女大奶头黄色视频| 99国产综合亚洲精品| 秋霞伦理黄片| 国产片特级美女逼逼视频| 人人澡人人妻人| av又黄又爽大尺度在线免费看| 少妇人妻久久综合中文| 国产成人精品无人区| 国产成人精品无人区| 久久久国产精品麻豆| 久久久国产精品麻豆| 久久婷婷青草| 国产乱来视频区| 春色校园在线视频观看| 国产精品久久久久久av不卡| 国产成人av激情在线播放 | 最黄视频免费看| 日本与韩国留学比较| 一级毛片黄色毛片免费观看视频| 亚洲国产精品999| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 亚洲国产日韩一区二区| 日本黄色日本黄色录像| 赤兔流量卡办理| 欧美97在线视频| 久久久久久久亚洲中文字幕| 十八禁高潮呻吟视频| 亚洲av综合色区一区| 一边摸一边做爽爽视频免费| 美女国产高潮福利片在线看| 一区二区三区乱码不卡18| 高清不卡的av网站| 午夜免费鲁丝| 国精品久久久久久国模美| 亚洲成人手机| 能在线免费看毛片的网站| 久久久久精品性色| 欧美3d第一页| 久久久久久久国产电影| 久久精品久久久久久久性| 免费看av在线观看网站| 久久精品久久久久久久性| 一边摸一边做爽爽视频免费| 国产成人aa在线观看| 亚洲精品第二区| 免费黄色在线免费观看| 日韩一区二区三区影片| 特大巨黑吊av在线直播| 精品人妻一区二区三区麻豆| 极品少妇高潮喷水抽搐| 在线观看免费视频网站a站| 成年人免费黄色播放视频| 中文字幕人妻丝袜制服| 夜夜看夜夜爽夜夜摸| 蜜桃在线观看..| 亚洲三级黄色毛片| 又黄又爽又刺激的免费视频.| 日本91视频免费播放| 日日啪夜夜爽| 午夜免费男女啪啪视频观看| 伦精品一区二区三区| 久久狼人影院| 亚洲欧美色中文字幕在线| 新久久久久国产一级毛片| 久久久久久久久久久丰满| 午夜福利视频在线观看免费| 国产精品麻豆人妻色哟哟久久| 热99久久久久精品小说推荐| 色5月婷婷丁香| 国产欧美亚洲国产| 亚洲国产精品一区二区三区在线| 街头女战士在线观看网站| 视频在线观看一区二区三区| 亚洲精品,欧美精品| 成人手机av| 国产一区二区三区av在线| 欧美日韩成人在线一区二区| 亚洲av国产av综合av卡| 亚洲av中文av极速乱| 最新的欧美精品一区二区| 国产午夜精品久久久久久一区二区三区| 久久99蜜桃精品久久| 国产精品成人在线| 亚洲av中文av极速乱| 久久99精品国语久久久| av播播在线观看一区| 91aial.com中文字幕在线观看| 97超视频在线观看视频| 欧美国产精品一级二级三级| 亚洲精品久久午夜乱码| 久久久久网色| 青春草国产在线视频| 纯流量卡能插随身wifi吗| 久久久久久久亚洲中文字幕| 人妻 亚洲 视频| 亚洲久久久国产精品| www.av在线官网国产| 老司机影院毛片| 狠狠婷婷综合久久久久久88av| 亚洲欧美一区二区三区黑人 | 丰满迷人的少妇在线观看| 久久国产亚洲av麻豆专区| 一级毛片我不卡| 久久久久精品性色| 日韩一区二区三区影片| 老司机影院成人| 日韩亚洲欧美综合| 国产亚洲精品久久久com| 久久精品夜色国产| 女的被弄到高潮叫床怎么办| av在线app专区| 日韩欧美一区视频在线观看| 国产精品.久久久| 人妻夜夜爽99麻豆av| 亚洲精品,欧美精品| 国产精品人妻久久久影院| 久久韩国三级中文字幕| 国产精品久久久久成人av| 亚洲国产精品成人久久小说| 啦啦啦啦在线视频资源| 99热这里只有精品一区| 亚洲av成人精品一区久久| a级片在线免费高清观看视频| 色94色欧美一区二区| 日韩av不卡免费在线播放| 成人黄色视频免费在线看| 亚洲av.av天堂| 人人妻人人添人人爽欧美一区卜| 国精品久久久久久国模美| a级毛片黄视频| 好男人视频免费观看在线| 黑人猛操日本美女一级片| 亚洲色图综合在线观看| 一本一本综合久久| 国产精品久久久久久精品电影小说| 一边摸一边做爽爽视频免费| 亚洲第一av免费看| .国产精品久久| 国产精品偷伦视频观看了| av又黄又爽大尺度在线免费看| 免费久久久久久久精品成人欧美视频 | 青春草亚洲视频在线观看| 蜜桃国产av成人99| 亚洲精品成人av观看孕妇| 纯流量卡能插随身wifi吗| 国产色爽女视频免费观看| 人妻 亚洲 视频| 最近中文字幕2019免费版| 美女大奶头黄色视频| 国产一区亚洲一区在线观看| 观看美女的网站| 黄色毛片三级朝国网站| 国产在视频线精品| 亚洲四区av| 青春草国产在线视频| 欧美成人精品欧美一级黄| 中文字幕久久专区| 亚洲精华国产精华液的使用体验| 亚洲国产精品一区二区三区在线| 女人精品久久久久毛片| 精品久久久久久久久av| 99九九在线精品视频| 亚洲精品成人av观看孕妇| av在线app专区| 免费高清在线观看日韩| 人妻夜夜爽99麻豆av| 一区在线观看完整版| 亚洲综合色网址| 一级二级三级毛片免费看| 一区二区三区免费毛片| 日韩 亚洲 欧美在线| 国产高清三级在线| 五月开心婷婷网| 日本午夜av视频| 在线观看三级黄色| 黄片播放在线免费| 精品一区二区三区视频在线| 亚洲少妇的诱惑av| 欧美激情国产日韩精品一区| 日本av免费视频播放| 国产一区二区在线观看av| 免费久久久久久久精品成人欧美视频 | 中文字幕亚洲精品专区| 国产又色又爽无遮挡免| 美女大奶头黄色视频| 国产欧美另类精品又又久久亚洲欧美| 色视频在线一区二区三区| 国产精品久久久久久精品电影小说| 高清毛片免费看| 国产又色又爽无遮挡免| 色网站视频免费| 男人操女人黄网站| 国产精品三级大全| 成人无遮挡网站| 欧美丝袜亚洲另类| 久久狼人影院| 亚洲成色77777| 超色免费av| 亚洲精品乱久久久久久| 人人妻人人添人人爽欧美一区卜| 一本一本综合久久| 久久久久久久大尺度免费视频| 男人添女人高潮全过程视频| 亚洲中文av在线| 一级,二级,三级黄色视频| 国产一区亚洲一区在线观看| 午夜福利影视在线免费观看| 中文字幕制服av| 不卡视频在线观看欧美| 免费高清在线观看视频在线观看| 国产男女内射视频| 国产高清有码在线观看视频| 国产精品国产三级国产专区5o| 黄色欧美视频在线观看| 欧美三级亚洲精品| 日韩精品有码人妻一区| 欧美另类一区| 亚洲综合精品二区| 亚洲精品乱码久久久v下载方式| av有码第一页| 视频在线观看一区二区三区| 97超视频在线观看视频| 国产av国产精品国产| 三级国产精品片| 另类亚洲欧美激情| 我的老师免费观看完整版| 亚洲av.av天堂| 一区二区三区四区激情视频| 亚洲国产精品国产精品| 久久久久久久久久人人人人人人| 国产黄频视频在线观看| 有码 亚洲区| 亚洲精品日本国产第一区| 99久久精品一区二区三区| 制服丝袜香蕉在线| 各种免费的搞黄视频| 久久午夜福利片| 97在线人人人人妻| 久久久久国产网址| 成人手机av| 国产成人av激情在线播放 | 国产精品蜜桃在线观看| 午夜福利,免费看| 亚洲美女黄色视频免费看| 制服人妻中文乱码| 91久久精品国产一区二区成人| 在线 av 中文字幕| 国产一区有黄有色的免费视频| 亚洲精品美女久久av网站| 日产精品乱码卡一卡2卡三| 亚洲av福利一区| 亚洲成人手机| av线在线观看网站| 日韩人妻高清精品专区| 国产高清国产精品国产三级| 久久精品国产a三级三级三级| 亚洲av日韩在线播放| 自拍欧美九色日韩亚洲蝌蚪91| 久久97久久精品| 亚洲国产成人一精品久久久| 色94色欧美一区二区| 精品久久久久久电影网| 好男人视频免费观看在线| 亚洲欧美色中文字幕在线| 欧美精品一区二区大全| 免费久久久久久久精品成人欧美视频 | 午夜免费观看性视频| 精品国产露脸久久av麻豆| 极品人妻少妇av视频| 亚洲成人手机| 久久精品国产鲁丝片午夜精品| 人体艺术视频欧美日本| 一级a做视频免费观看| 人妻系列 视频| 日韩熟女老妇一区二区性免费视频| 国产国语露脸激情在线看| 国精品久久久久久国模美| 亚洲精品乱码久久久v下载方式| 久久精品人人爽人人爽视色| 亚洲综合色惰| 国内精品宾馆在线| 日本欧美视频一区| 男女啪啪激烈高潮av片| 久久久久久伊人网av| videossex国产| 精品视频人人做人人爽| 久久久亚洲精品成人影院| 国产亚洲最大av| 国产精品久久久久久精品古装| 国产男女超爽视频在线观看| 亚洲精品国产色婷婷电影| 亚洲av在线观看美女高潮| 一级,二级,三级黄色视频| 国产亚洲最大av| 超碰97精品在线观看| 久久鲁丝午夜福利片| 日本wwww免费看| 一级毛片我不卡| 亚洲精品一二三| 日日撸夜夜添| 久久精品夜色国产| 国产伦理片在线播放av一区| 汤姆久久久久久久影院中文字幕| 视频中文字幕在线观看| 建设人人有责人人尽责人人享有的| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | av免费观看日本| 国产高清有码在线观看视频| 99热这里只有精品一区| 欧美最新免费一区二区三区| 亚洲一级一片aⅴ在线观看| 一本大道久久a久久精品| 在线观看一区二区三区激情| 大话2 男鬼变身卡| 2018国产大陆天天弄谢| 亚洲国产色片| 亚洲人成网站在线观看播放| 国产精品偷伦视频观看了| av国产精品久久久久影院| 内地一区二区视频在线| 99国产综合亚洲精品| 国产成人精品婷婷| 性色av一级| 波野结衣二区三区在线| 成人国产av品久久久| 97超碰精品成人国产| 国产精品一国产av| 大片免费播放器 马上看| 亚洲国产毛片av蜜桃av| 欧美精品一区二区大全| 国产亚洲最大av| 成人二区视频| 老司机影院毛片| 男女边吃奶边做爰视频| 99热这里只有是精品在线观看| www.av在线官网国产| 亚洲综合色网址| 亚洲国产欧美日韩在线播放| 亚洲欧美一区二区三区黑人 | 91精品三级在线观看| 十分钟在线观看高清视频www| 亚洲欧美一区二区三区国产| 国产高清不卡午夜福利| 午夜久久久在线观看| 内地一区二区视频在线| 22中文网久久字幕| 啦啦啦中文免费视频观看日本| 中文字幕久久专区| 美女主播在线视频| 国产精品久久久久久精品电影小说| 欧美老熟妇乱子伦牲交| 99久久中文字幕三级久久日本| 人妻系列 视频| 国产成人精品在线电影| 成人18禁高潮啪啪吃奶动态图 | 亚洲精华国产精华液的使用体验| 搡女人真爽免费视频火全软件| 亚洲国产色片| 免费看光身美女| 又粗又硬又长又爽又黄的视频| 成人二区视频| 国产成人精品久久久久久| 国产综合精华液| 久久久国产欧美日韩av| av网站免费在线观看视频| 观看美女的网站| videos熟女内射| 蜜臀久久99精品久久宅男| 九色成人免费人妻av| 国产视频内射| 性高湖久久久久久久久免费观看| 超色免费av| av电影中文网址| 97超碰精品成人国产| 国产精品秋霞免费鲁丝片| 国产视频内射| 成人毛片60女人毛片免费| xxxhd国产人妻xxx| 久久精品熟女亚洲av麻豆精品| 少妇丰满av| 国产 精品1| 亚洲人成77777在线视频| 国产片内射在线| 亚洲人成网站在线播| 三级国产精品欧美在线观看| 午夜福利视频精品| 国产成人freesex在线| 欧美精品国产亚洲| 一个人免费看片子| 丝瓜视频免费看黄片| 欧美三级亚洲精品| 美女xxoo啪啪120秒动态图| 亚洲人成网站在线播| 永久免费av网站大全| av在线播放精品| 久久狼人影院| 美女主播在线视频| 久久久久网色| 欧美日韩综合久久久久久| 久久综合国产亚洲精品| 成年女人在线观看亚洲视频| 亚洲av成人精品一二三区| 久久精品夜色国产| 国产精品久久久久久久电影| 一二三四中文在线观看免费高清| 两个人免费观看高清视频| 一区二区av电影网| 青春草国产在线视频| 中文字幕人妻丝袜制服| 日本爱情动作片www.在线观看| 日韩欧美一区视频在线观看| 国产成人91sexporn| 久久99热6这里只有精品| 国产av码专区亚洲av| 91精品国产国语对白视频| 国产av一区二区精品久久| 一区二区av电影网| 国产精品99久久99久久久不卡 | 赤兔流量卡办理| 国产一区二区在线观看av| 国产视频内射| 男女免费视频国产| 最后的刺客免费高清国语| kizo精华| 2018国产大陆天天弄谢| 自线自在国产av| 蜜臀久久99精品久久宅男| 中文字幕精品免费在线观看视频 | 麻豆成人av视频| 一级毛片 在线播放| 久久精品久久久久久久性| 大话2 男鬼变身卡| 能在线免费看毛片的网站| a级片在线免费高清观看视频| 国产精品免费大片| 欧美成人午夜免费资源| 91国产中文字幕| 人妻少妇偷人精品九色| 国产av精品麻豆| 色视频在线一区二区三区| 亚洲不卡免费看| 国产精品成人在线| 国产精品99久久99久久久不卡 | 麻豆成人av视频| 97在线视频观看| 久久精品国产鲁丝片午夜精品| 超色免费av| 少妇 在线观看| 久久女婷五月综合色啪小说| 国产精品一区二区在线不卡| 国产淫语在线视频| 国产精品 国内视频| 日本欧美国产在线视频| 免费大片黄手机在线观看| 激情五月婷婷亚洲| 青春草视频在线免费观看| 桃花免费在线播放| 国产深夜福利视频在线观看| 久久久久精品久久久久真实原创| 丰满饥渴人妻一区二区三| 亚洲精品aⅴ在线观看| 高清欧美精品videossex| 伦理电影大哥的女人| 天美传媒精品一区二区| 国产不卡av网站在线观看| 高清av免费在线| 少妇精品久久久久久久| 日本vs欧美在线观看视频| 三上悠亚av全集在线观看| 黄片无遮挡物在线观看| 国产亚洲精品久久久com| 久久女婷五月综合色啪小说| 黄色怎么调成土黄色| 久久久欧美国产精品| 又粗又硬又长又爽又黄的视频| videos熟女内射| av免费在线看不卡| 亚洲美女搞黄在线观看| 久久久久精品久久久久真实原创| 欧美 日韩 精品 国产| 免费人妻精品一区二区三区视频| 最近2019中文字幕mv第一页| 亚洲国产精品一区二区三区在线| 性色avwww在线观看| 亚洲成色77777| av卡一久久| 看免费成人av毛片| 午夜免费鲁丝| 国产毛片在线视频| 久久午夜综合久久蜜桃| 99久久精品一区二区三区| 777米奇影视久久| 午夜影院在线不卡| 欧美亚洲 丝袜 人妻 在线| 99九九在线精品视频| 亚洲经典国产精华液单| 亚洲内射少妇av| 美女视频免费永久观看网站| 国产精品成人在线| a级毛色黄片| 18在线观看网站| 三级国产精品片| 2022亚洲国产成人精品| 51国产日韩欧美| 欧美精品人与动牲交sv欧美| a级毛片黄视频| 边亲边吃奶的免费视频| 曰老女人黄片| 人成视频在线观看免费观看| 在线 av 中文字幕| 国产免费福利视频在线观看| 大香蕉97超碰在线| 亚洲欧美精品自产自拍| 99久久精品国产国产毛片| 满18在线观看网站| 丝袜脚勾引网站| 亚洲精品久久久久久婷婷小说| 蜜桃国产av成人99| 国产黄色免费在线视频| 青春草国产在线视频| 免费高清在线观看视频在线观看| 少妇熟女欧美另类| 国产精品一国产av| 在线亚洲精品国产二区图片欧美 | 黄色视频在线播放观看不卡| 欧美精品人与动牲交sv欧美| 国产黄色视频一区二区在线观看| 母亲3免费完整高清在线观看 | 日本黄色日本黄色录像| 丝袜美足系列| 肉色欧美久久久久久久蜜桃| 成年人午夜在线观看视频| 一级,二级,三级黄色视频| 国产亚洲一区二区精品| 日本色播在线视频| 丝袜在线中文字幕| 亚洲av成人精品一二三区| 亚洲图色成人| 一本色道久久久久久精品综合| 亚洲久久久国产精品| 在线观看免费高清a一片| 免费av不卡在线播放| 日韩免费高清中文字幕av| 国产免费福利视频在线观看|