999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

長江口南北港水下地形演變數值模擬研究

2023-04-30 17:18:45郭佳李昌文楊忠勇孫詩為欒華龍
人民長江 2023年13期
關鍵詞:區域模型

郭佳 李昌文 楊忠勇 孫詩為 欒華龍

摘要:

南北港是長江口“三級分叉、四口入海”格局中的二級河道,也是連接南支河道和北槽深水航道的重要通航水路,對其沖淤演變特征及未來演變趨勢開展數值研究具有重要意義。首先構建了包含整個長江口的水沙動力學數值模型,在實測資料驗證下使用地貌加速因子對未來60 a間南北港區域的河道演變趨勢進行了模擬分析,同時結合南北港2006~2016年間的實測水深數據開展對比分析。結果表明:① 在新的來水來沙條件以及本地青草沙工程影響下,南北港水下地形經歷了顯著的沖淤演變過程,實測數據和數值模擬結果均顯示出“灘淤槽沖”的基本演變格局,即河道沖淤深度與水深基本呈線性關系;② 沖淤演變的轉換水深約為-4~-6 m,未來60 a內深槽區域的最大沖刷深度可達-12 m;③ 基于河口泥沙的空間運動滯后效應,分析此種河口“灘淤槽沖”模式實際上是由岸灘和深槽之間水質點和泥沙顆粒輸運的相位差所致。

關 鍵 詞:

長江口南北港; 沖淤演變; 灘淤槽沖; 空間滯后效應; 數值模擬

中圖法分類號: TV82

文獻標志碼: A

DOI:10.16232/j.cnki.1001-4179.2023.S2.033

0 引 言

長江是世界第三長河,長江流域巨大的輸沙量促成了長江口豐富的灘涂資源,不過近年來在全球氣候變化作用以及上游河道水利工程的影響下,長江口水下地貌也發生了復雜變化。長江口復雜的來水條件和來沙條件,以及長江口面臨的徑流、潮流相互作用,是長江口水下地貌發生復雜演變的重要原因[1-2]。陳吉余[3]、楊世倫[4]、戴志軍[5]等通過系統研究,將長江口2000年以來的河道發育和演變模式總結為“南岸邊灘推展,北岸沙島并岸,河口束狹,河道成形,河槽加深”,其研究成果為后來河口演變方面的研究和圍海造陸等灘涂資源綜合開發與治理提供了指導。

大通站統計資料顯示,1960年以來長江口年均徑流量約8 906億m3/a,雖然時有洪水發生,但來水量的年際變化特征并不明顯,不過徑流量在年內有明顯的季節性變化特征,汛期主要集中在每年的5~10月,徑流量占比約70.3%,每年11月至次年4月為非汛期。長江流域的來沙豐富,大都以細顆粒泥沙為主,泥沙粒徑一般在0.04 mm以下,長江口口門區域由于水流突然展寬,加之外海潮流的作用,因此泥沙在此沉降,形成攔門沙區域,同時也形成了大量的潮汐灘涂濕地。自20世紀70年代以來,由于長江干支流開始建壩,同時長江流域也實施了大量的退耕還林政策,長江輸沙量逐漸減少。統計表明,1960~1968年間,長江年均輸沙量為5.22億 t,1985年流域輸沙量略微降低至4.49億 t,到2000年時流域輸沙量降低至3.49億 t。2004年后,隨著三峽水庫蓄水,輸沙量發生銳減,年輸沙量降低至1.47億 t,不足1960年左右的30%[6-7]。

長江口來水和來沙量發生巨變后,河口三角洲水下地貌形態及其穩定性必然發生顯著變化。大量研究顯示,三峽工程的修建是長江口來沙量發生劇減的轉折點[8],此外長江各個支流上修建大壩以及國家的退耕還林政策實施也是導致來沙量降低的重要原因[9-10]。不過Dai等[11]研究發現,在長江上游來沙量劇減的大背景下,長江口的水下地貌特征自2013年起已經逐漸轉為沖淤平衡狀態。長江口北港作為連接北槽深水航道的重要汊道,對長江口的航運暢通[12-13]及綜合治理[14-15]起著重要作用。因此針對新形勢下長江口北港地貌演變開展研究具有重要意義。

1 研究區域概況

長江口地形態勢如圖1所示,長江口自徐六涇開始,下游河段被崇明島分隔成南支和北支兩個汊道,南支河道又進一步被長興島和橫沙島分隔成南港和北港兩個二級汊道,繼續往下游,南港河道再一次被九段沙分割為南槽和北槽兩個三級汊道,最后長江口以“北支-北港-北槽-南槽”四口入海的形式與大海交接。由于北槽河道有重要工程深水航道,南港和北港是連接深水航道與上游南支河道的重要紐帶,同時南港河道中暗沙較多,河道演變頻繁,發育有青草沙、中央沙等大型沙體,因此南港和北港的地貌演變十分重要。南北港的河勢演變主要受上游南支的與北支的分流和分沙影響是最大的,歷史上北支為主要出海口時,南港和北港的面積都非常小。1954年長江發生特大洪水后,長江南支的主汊經過崇明島水道進入南港,在此期間大量泥沙淤積在南港河道形成瑞豐沙等,從此南港由單一河槽變為復式河槽。長江口主要有9個潮位站點,包括青龍港(QLG)、連興港(LX)、堡鎮(BZ)、吳淞(WS)、橫沙(HS)、中浚(ZJ)、九段東(JDD)、雞骨礁(JGJ)和小洋山(XYS),以及6個潮流站點,包括南支(SB)、南港(SC)、南槽(SP)、長興(CX)、北港(NC)和北槽(NP),所有站點的具體位置如圖1所示。

2 數值模型構建及驗證

為分析長江口南北港河段地貌演變過程及影響因素,本文采用數值模擬進行研究。MIKE21是一種專業模擬河道水流、泥沙和地貌演變等的數值模型,模型所采用的有限體積非結構網格能夠很好地保證通量守恒,而且非結構三角形網格可以很好地擬合島嶼等復雜岸線,提高模擬的精度。因此本次研究中選擇MIKE21 FM開展數值模擬計算。

2.1 模型區域、網格劃分及邊界條件

本次研究中構建的數值模型范圍如圖1所示,包括整個長江口潮流界以下范圍,杭州灣和部分東海近海岸的區域。長江口上游設在潮流界江陰以上約100 km處,保證雙向潮流不受上游邊界的影響。杭州灣上游入口設置在約120.2 °E處,以實際的水位波動給定。外海開邊界上,北部海洋開邊界設置在約32.3 °N處,南部海洋開邊界設置在約29.2 °N處,東部海洋開邊界設置在約124 °E處,各海洋邊界均以實際水位波動給定。整個模型范圍內,共劃分約72 100個節點和141 680個網格(圖1左下角顯示了南港和北港部分網格)。長江上游開邊界區域,網格分辨率約600 m;外海開邊界區域,網格分辨率放大至約5 000 m;在本次研究重點關心的區域(即長江口范圍內),網格分辨率約300 m。

模型中長江口范圍內的水深地形數據均來源于最新的海圖測量數據。長江口外的東海區域水深數據來源于NOAA網上下載的地形數據(http:∥www.ngdc.noaa.gov)。河流邊界的徑流量和輸沙量采用大通水文站的實測數據。由于東海前進潮波主要從東南方向傳入長江口,故將東邊界和南邊界作為潮波入射邊界,設置8 個主要分潮(M2、S2、N2、K2、K1、O1、P1、Q1),各個分潮的調和常數由中國海洋局部模型計算成果提供(http:Volkov.oce.orst.edu/tides),并采用《渤海黃海東海海洋圖集(水文)》的數據進行矯正,從而計算相應時刻的水位。模型范圍內,各區域的初始懸沙濃度采用2013年泥沙實測數據。

2.2 模型水位驗證

針對模型計算的水位,主要選取長江口9個潮位站點(具體位置見圖1)的調和常數(振幅An和位相Gn)以及實測水位進行驗證。將模型計算的9個潮位站點的水位進行調和分析計算,與各個分潮潮汐表上的潮位分析結果進行對比。結果表明,部分站點模型計算出現波峰的時間較實測略有延遲,但計算所得的水位與實測水位的振幅和周期基本一致,吻合度較高。綜上所述,本文采用的模型計算結果具有較高的可信度。

2.3 模型流速流向驗證

根據前文構建的數值模型,分別針對圖1中所示6個站點上表層水體的流速和流向進行對比,根據監測資料,模型驗證時間為2018年5月5日02:00至5月7日03:00,該時段為天文大潮期間,長江口內潮流流速垂線平均值最大可達2.0 m/s。驗證結果表明:各站點流速計算值與實測值差異較小,僅在部分站點的高流速期間計算值偏大,而實測值偏小,但這并不影響整體結果。從流向計算結果上看,計算結果和實測結果基本吻合,不過在漲落潮流轉向期間,流向模擬結果誤差較大,主要原因是該期間流速小,流場較紊亂。整體來講,本次構建并調試過后的數值模型模擬精度較高,能夠滿足研究所需。

3 長江口南北港沖淤現狀及趨勢分析

3.1 水深特征

2006~2016年,長江口南北港河槽以及周邊河道水下地貌特征如圖2所示。由于2006年和2009年獲取的實測數據分辨率較低(約100 m),因此得到的水下地形圖不如2011年和2016年的平滑。總體來說,各年間水下地形模式基本一致,這10 a來在水深上并未發生顯著變化。南北港河道主流區域水深較深,最深處可達20 m左右,且這一深度隨著時間演進而逐漸增大。南北港河道在2009年發生較為顯著的演變特征,2006年時水深約為15 m左右,2009年由于處于青草沙水庫施工期間,水體含沙量本應相對較高,泥沙落淤后水深較淺,上口段平均水深不足10 m。不過2011年后由于北港上口段束窄后流速變大,強勁的水動力沖刷河槽,河槽加深,最大水深又重新恢復至16 m,2011年后北港河道水深逐漸趨于穩定,從側面論證了工程影響下,局部地形發生劇烈改變但又迅速調整恢復的過程。

3.2 沖淤特征

根據2006年、2011年和2016年長江口南北港水下地貌實測數據,繪制研究區域水下地貌沖淤演變趨勢如圖3所示。河道地形特征和沖淤演變趨勢,與斷面懸沙分布特征和輸運結構密切相關,圖中正值表示淤積,負值表示沖刷。近年來北港為微彎型河勢,其河道中央深槽主要分布在中偏南側。在此微彎河道橫斷面上形成逆時針環流結構,近底層水體從南側向北側輸運,導致泥沙大量分布在河槽北岸。從輸沙結構特征上來說,北港主槽河道北側大部分水體總輸沙沖向大海、南側小部分輸沙方向與漲潮流方向一致,上溯泥沙在北港上段遇落潮流而落淤,將進一步沖刷河槽北岸,導致微彎趨勢逐漸增強。2009年南北港分汊口控制工程以及青草沙水庫等大型工程修建后,北港上口通道演變得以控制,河道束窄,河道的彎曲程度呈增大的趨勢。此外,橫沙通道的河道走向致使漲潮流直指崇明島南岸,也是研究斷面南側深槽形成原因之一。可見圖3(a)中北港主槽橫斷面泥沙分布特征和輸沙結構特征與河勢演變趨勢是基本一致的。圖3(b)所示的近年北港沖淤趨勢主要呈“主槽沖刷,兩岸淤積”的特點,一方面可歸因于北港微彎河勢進一步加強,主槽區域進一步刷深;另一方面與泥沙自身運動特征中的空間滯后效應相關。河口地貌在適應上游水沙條件及大型工程影響下的自適應調整中,動力較強的深槽區域泥沙易于侵蝕難以沉降,而動力較弱的淺灘區域泥沙易于沉降難以侵蝕[16],導致“灘淤槽沖”的現象出現,這種現象在其他入海汊道如北港同樣存在[1]。

3.3 未來60 a南北港河道沖淤演變趨勢

根據數值模型中的地貌加速因子,計算未來60 a間北港河道沖淤趨勢分布,模型加速時間比尺為1∶365,即模型計算1 d相當于1 a,結果如圖4所示。圖中暖色調(黃色和紅色)表示淤積區域,冷色調(青色和綠色)表示沖刷區域,白色等值線表示沖淤值為0的區域。圖4(a)中顯示了未來10 a長江口的沖淤情況模擬結果。從圖中可以看出,未來10 a間,沖淤格局的分布已經呈現出較為明顯的趨勢。沖刷區域集中在各個過水的通道,例如南支和南港的南部深槽區域、北支下段口門區域、北港北沙以南的區域、九段沙南部的南槽中部深槽區域和北部的深水航道下段區域,但沖刷幅度較小,最大沖刷值僅約-2 m左右。淤積區域集中在各個灘涂區域,例如南支北部的東風西沙區域、北支上游河段、崇明東灘區域、橫沙島以東的區域、九段沙區域以及南匯邊灘區域,淤積幅度也較小,最大淤積值僅2 m左右。整個河口的沖淤格局基本呈現“灘淤槽沖”的趨勢,在不考慮人類工程和其他外來因素的影響條件下,沖淤幅度較低,凈沖淤值僅約0.2 m/a。

20 a后,長江口的沖淤格局和沖淤區域與10 a內基本一致,但沖淤的幅度略有加強,最大的沖刷和淤積幅度值都在2 m以上(見圖4(b))。從20 a后長江口的水下地貌與初始地貌對比分析來看,-5 m等深線似乎進一步向陸地方向移動,其包圍的區域逐漸在減小,-2 m等深線與初始水下地貌時相比幾乎維持不動。然而,0 m等深線在逐漸向外海方向移動,即0 m等深線以上的面積有擴大的趨勢,這也正是“灘淤槽沖”的結果。

圖4(c)~4(f)分別顯示了未來30 a至未來60 a長江口的沖淤圖和水下地貌。可以看出,各個10 a間,長江口中“灘淤槽沖”的格局基本沒有改變。各個淺灘區域總是發生不同程度的淤積,各個深槽區域總是發生不同程度的沖刷。隨著模擬時間的加長,沖刷和淤積的幅度在逐漸變大,但沖刷幅度增長較快,由于淺灘區域水深本來就小,可供淤積的空間較小,因此淤積的幅度增長較慢。到60 a后,最大的沖刷幅度值可達-6~-8 m,最大的淤積幅度可達4~6 m。從60 a后各等深線的演變情況來看,-5 m等深線顯著地往內陸移動,包圍區域的面積顯著減小,0 m等深線則顯著地往外海方向移動,包圍的面積顯著增加,-2 m等深線則變化不大,其包圍區域的面積略有增加。

總體而言,沿著主泓線方向河槽有明顯的加深趨勢,崇明島南沿以及青草沙水庫北側略有淤積,北港河道基本趨于穩定。未來為維持北港兩側岸線穩定,建議需重點加固進口段北側岸線以及橫沙島北沿的岸灘。若將-5 m等深線以內的區域都視為潛在可圍墾區域的話,在未施加任何人工措施的情況下,未來幾十年內,其面積是逐漸減小的,不過減小的幅度不大。但是0 m等深線以上的區域面積在未來幾十年內會逐漸擴大,表明易于圍墾區域的面積還在逐漸增加。

4 討 論

從前文中針對南北港沖淤現狀特征及趨勢的分析可知,在長江上游大量建壩及實施退耕還林等政策影響下,長江來水來沙量發生顯著變化,南北港河道基本呈現“深槽區域沖刷、淺灘區域淤積”的特征。為此,本文根據圖3~4中沖淤演變數據,統計了不同水深區域的沖淤特征,結果如圖5所示。不管是實測數據分析結果還是未來預測結果,均顯示出顯著的“灘淤槽沖”演變趨勢。在水深20 m區域的河道,2006~2016年平均沖刷深度約-7.3 m,未來60 a可達-11.2 m。而當水深超過-3 m后,實測值和預測值均顯示出淤積特征,不過淤積幅度較小,均不足2 m。灘淤槽沖的轉換深度約在-7~-4 m之間。

河口的這種“灘淤槽沖”現象從理論上可用“空間沉降滯后”效應進行解釋。圖6描述了經典空間沉降滯后效應下水質點和泥沙顆粒運動示意。圖6是根據Nichols等[17]的成果改繪而來,圖中實線表示水質點運動軌跡,虛線表示泥沙顆粒運動軌跡。河口泥沙空間滯后效應最早由Postma[18]在Wadden海域的灘槽水沙輸運研究中注意到,他發現由于平均潮流速在“灘→槽”方向遞增,使得漲潮過程中的向岸輸沙量大于落潮過程中的離岸輸沙量。該模式經Groen[19]進一步總結并修正。漲潮過程中,泥沙顆粒在水質點A的作用下于位置①起動,并向岸邊輸運,到位置②時逐漸沉降,落淤于位置③。落潮過程中,水質點A不能再起動沉降于位置③的泥沙顆粒,而由靠岸的水質點B起動,最后在位置④開始沉降,并落淤于⑤點,如此一個潮周期后,泥沙顆粒將發生向岸的凈位移(①→⑤),即為泥沙的空間滯后效應所致。考慮到非均勻沙的影響,泥沙輸運的空間滯后效應實際上是潮流環境中水沙質點運動的歐拉-拉格朗日不對稱性導致的泥沙在空間上的輸移滯后過程,也就是泥沙質量守恒方程中泥沙顆粒和水流質點相互作用的非線性對流項所貢獻的含沙量和輸沙量[20]。

長江口南北港這種灘淤槽沖現象在之前的研究中也有描述[2]。通過近30 a來長江口多個河槽(包括北支、北港、南匯等)的岸灘剖面形狀及淺地層沉積物演變特征分析發現,岸灘剖面呈現坡腳刷深、坡頂淤積的現象。Luan等[21]通過數值模型預測了長江口未來20 a的水下地貌演變情況,其預測結果和長江口多個河槽的實測資料分析結果[22]均顯示出“灘淤槽沖”的演變格局。

5 結 論

長江口南北港上承南支河道,下接長江口深水航道,是長江口重要的通航水路。本文根據2006~2016年間長江口南北港區域的實測水深數據分析了10 a間該區域的沖淤演變特征,同時基于數值模型模擬,預測了未來60 a間南北港的沖淤演變趨勢,主要結論如下。

(1) 長江口南北港河道在2009~2016年間經歷了較為顯著的沖淤演變過程,特別是2009年青草沙水庫建設期間,長興島頭部區域淤積深度近6 m,不過隨著工程建設的結束,深槽重新形成。

(2) 數值模擬結果顯示,未來60 a間長江口基本維持“灘淤槽沖”的演變趨勢,-5 m等深線顯著地往內陸移動,包圍區域的面積顯著減小;0 m等深線則顯著地往外海方向移動,包圍的面積顯著增加;-2 m等深線則變化不大,其包圍區域的面積略有擴大。

(3) 實測數據和數值模擬結果的統計結果表明,長江口南北港沖淤深度與水深基本呈線性關系,即水深越大,沖刷越強,初步推斷這種“灘淤槽沖”模式是長江口泥沙輸運過程中的空間滯后效應所致。

參考文獻:

[1] 郭興杰,程和琴,楊忠勇,等.長江口北港河勢演變及趨勢分析[J].泥沙研究,2016(5):33-39.

[2] 計娜,程和琴,楊忠勇,等.近30年來長江口岸灘沉積物與地貌演變特征[J].地理學報,2013,68(7):937-946.

[3] 陳吉余,惲才興,徐海根.兩千年來長江河口的發育模式[J].海洋學報,1979,1(1):103-111.

[4] 楊世倫,吳秋原,黃遠光.近40年崇明島周圍灘涂濕地的變化及未來趨勢展望[J].上海國土資源,2019,40(1):68-71.

[5] 戴志軍.海岸與近海環境變化過程[J].熱帶海洋學報,2022,41(4):3-4.

[6] 萬遠揚,吳華林.徑流量變化對長江口北槽最大渾濁帶影響分析[J].水利水運工程學報,2021(5):1-7.

[7] 左書華,楊春松,付桂,等.長江口入海水沙通量變化及其影響分析[J].海洋地質前沿,2022,38(11):56-64.

[8] GAO A,YANG S L,LI G,et al.Long-term morphological evolution of a tidal island under influences of natural episodes and human activities,the Yangtze Estuary[J].Journal of Coastal Research,2010,26(1):123-131.

[9] YANG S L,LUO X X,STIJN T,et al.Role of delta-front erosion in sustaining salt marshes under sea-level rise and fluvial sediment decline[J].Limnology and Oceanography,2020,65(9):1990-2009.

[10] LUO X X,YANG S L,WANG R S,et al.New evidence of Yangtze delta recession after closing of the Three Gorges Dam[J].Scientific Reports,2017(7):41735.

[11] DAI Z J,LIU J T,WEI W,et al.Detection of the Three Gorges Dam influence on the Changjiang (Yangtze River)submerged delta[J].Scientific Reports,2014(4):6600.

[12] 吳帥虎,程和琴,鄭樹偉,等.近期長江口北槽河段河槽演變對人類活動的響應[J].長江流域資源與環境,2020,29(6):1401-1412.

[13] 應銘,季嵐,周海.長江口北槽12.5 m深水航道回淤的物理過程[J].水運工程,2017,536(11):77-85.

[14] 韓玉芳,竇希萍.長江口綜合治理歷程及思考[J].海洋工程,2020,38(4):11-18.

[15] 陳正兵,陳前海,侯衛國,等.新情勢下長江口綜合整治的思考[J].人民長江,2022,53(增 1):1-4.

[16] CHEN Y,LI Y,THOMPSON C,et al.Differential sediment trapping abilities of mangrove and saltmarsh vegetation in a subtropical estuary[J].Geomorphology,2018,318:270-282.

[17] NICHOLS M M,BIGGS R B.Estuaries:coastal sedimentary environments[M].New York:Springer,1985:77-186.

[18] POSTMA H.Hydrography of the Dutch Wadden Sea[J].Archives Néerlandaises De Zoologie,1954,10(4):405-511.

[19] GROEN P.On the residual transport of suspended matter by an alternating tidal current[J].Netherland Journal of Sea Research,1967,3(4):564-574.

[20] DE SWART H E,ZIMMERMAN J.Morphodynamics of tidal inlet systems[J].Annual Review of Fluid Mechanics,2009,41(1):203-229.

[21] LUAN H,DING P,WANG Z,et al.Process-based morphodynamic modeling of the Yangtze Estuary at a decadal timescale:controls on estuarine evolution and future trends[J].Geomorphology,2017,290:347-364.

[22] ZHANG X,LI J,ZHU W,et al.The self-regulation process and its mechanism of channels′ bed changes in the Changjiang (Yangtze)Estuary in China[J].Acta Oceanologica Sinica,2015,34(7):123-130.

(編輯:胡旭東)

猜你喜歡
區域模型
一半模型
永久基本農田集中區域“禁廢”
今日農業(2021年9期)2021-11-26 07:41:24
分割區域
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
關于四色猜想
分區域
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
基于嚴重區域的多PCC點暫降頻次估計
電測與儀表(2015年5期)2015-04-09 11:30:52
主站蜘蛛池模板: 996免费视频国产在线播放| 亚洲成人网在线播放| 国产精品亚洲片在线va| 欧美久久网| 99在线观看免费视频| 国产午夜福利亚洲第一| 成人午夜天| 伊人网址在线| 亚洲精品黄| 欧美无专区| 久久亚洲天堂| 无码人中文字幕| 亚洲精品国产首次亮相| 99在线观看视频免费| 香蕉视频在线精品| 午夜精品一区二区蜜桃| 91久久精品日日躁夜夜躁欧美| 一区二区欧美日韩高清免费| 朝桐光一区二区| 高清无码手机在线观看| 国产素人在线| 亚洲精品福利网站| 成人亚洲天堂| 亚洲av无码成人专区| 久久久久久午夜精品| 欧美日在线观看| 久草热视频在线| 欧美精品黑人粗大| 国产精品天干天干在线观看| 午夜天堂视频| 亚洲国产在一区二区三区| 在线免费不卡视频| 制服无码网站| 中文字幕2区| 日本中文字幕久久网站| 国产成人免费高清AⅤ| 国产精品无码一二三视频| 欧美全免费aaaaaa特黄在线| 成人久久18免费网站| 亚洲一级毛片免费观看| 色亚洲成人| 综合久久久久久久综合网| 国产精品高清国产三级囯产AV| 欧美亚洲另类在线观看| 2021天堂在线亚洲精品专区| 国产在线无码一区二区三区| 亚洲一区无码在线| 国产中文在线亚洲精品官网| 色国产视频| 国产成人综合网| 国产福利大秀91| 国产黄在线免费观看| 国产精品第一区| 国产乱视频网站| 亚洲综合婷婷激情| 99热这里只有精品5| 国产成人啪视频一区二区三区| 欧美a级在线| 情侣午夜国产在线一区无码| 午夜福利视频一区| 久久91精品牛牛| 麻豆精品久久久久久久99蜜桃| 欧美午夜久久| 东京热高清无码精品| 免费又爽又刺激高潮网址 | 最近最新中文字幕免费的一页| 97国产在线观看| 999福利激情视频| 狠狠ⅴ日韩v欧美v天堂| 2048国产精品原创综合在线| 亚洲综合色婷婷中文字幕| 亚洲午夜国产精品无卡| 一级爱做片免费观看久久| 亚洲 欧美 偷自乱 图片| 2021国产精品自产拍在线| 欧美日韩亚洲综合在线观看| 成·人免费午夜无码视频在线观看 | 国产丝袜啪啪| 国产69精品久久久久妇女| 好紧好深好大乳无码中文字幕| 97成人在线观看| www精品久久|