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

長江河口水流輸運時間研究

2014-04-18 12:01:08王亞何青沈健
海洋學報 2014年1期
關鍵詞:研究

王亞,何青*,沈健

(1.華東師范大學 河口海岸學國家重點實驗室,上海 200062;2.Virginia Institute of Marine Science,College of William and Mary,Gloucester Point,VA.,23062,USA)

1 引言

在水環境系統中,絕大多數生物、溶解營養鹽、污染物和懸浮顆粒均隨水體運動;在過去的幾十年,過量的營養鹽排放導致近岸海域、河口、湖泊的嚴重的環境問題。近30a來,人們普遍認識到水體物理過程和生物化學過程相互影響并控制污染物的歸宿,在近海水域生態系統問題中起著重要的作用[1-2]。因此,在研究河口環境問題之前,定量研究水體輸運過程至關重要。

多數海洋運動過程非常復雜,常規的現場觀測和數值模擬,多著眼于常規的狀態變量及其過程,難以衡量它們之間的內在機理,因此就需要一些特定的解讀方法,引入一些輔助性的變量來說明,水齡(water age)就是這一類的輔助變量。以水齡、滯留時間(residence time)和通過時間(transit time)為代表的時間尺度,提供了一個量化動力與系統功能之間聯系的可能性,這些變量不僅可以進一步理解模型計算結果,更重要的是在跨學科交叉研究特別是水環境管理方面是一個極其有益的方法[3-5]。

水交換一直是海洋學科中基礎的科學問題之一。目前,箱式、質點追蹤和保守物質的擴散模型在我國得到廣泛的應用,使得水體交換能力的研究得到長足發展。匡國瑞等[6]以高低潮鹽度變化給出了乳山灣一個潮周期內的水交換率。朱小兵等[7]依據全潮水文觀測資料,采用平均滯留時間研究海南島博鰲港水體交換能力。姬厚德等[8]利用實測水文資料采用納潮量的方法進行了筼筜湖海水平均半交換時間的研究。對于岸線和水深變化相對簡單區域,以上方法計算簡單,不失為有效的方法,但是對于岸線和水深變化復雜的區域,上述方法計算結果精度較差。水交換的另一重要方法保守性物質的擴散的方法常被用來計算河口、海灣的水交換,董禮先和蘇紀蘭[9]應用于象山港水交換的研究;Liu等[10]應用此方法研究膠州灣的平均滯留時間;王聰等[11]采用保守物質擴散方法研究平均風場作用下大亞灣水體交換能力。由于保守物質濃度存在復雜時空變化,在研究海灣河口的交換能力時,需要對時間和空間進行積分,所以此方法難以表達在河口海灣區域水流交換能力的時空差異,在近岸海域的水流交換能力研究方面存在局限性。

針對上述水交換方法的不足,科學家先后提出了一系列新的水交換的計算方法,其中最具有代表性并被廣泛采用的是“水齡”法。空間上某點的水齡定義為水體進入研究區域到該點所需要的時間;在Bolin與 Rohde[12]和 Takeoka[1-2]研究中,水體滯留時間定義為從該點到離開研究區域所需要的時間。近年來,這些概念被發展和應用在潟湖、河口、湖泊和深海的輸運時間的研究中。

Deleersnijder等[3]首先引入計算水齡的歐拉理論,并將其運用在英吉利海峽和北海南部平均水齡的空間分布。從此數值模擬在水齡和滯留時間的研究中得到廣泛的應用。Shen和 Hass[13]與Shen和Lin[4]分別利用水齡理論研究美國Chesapeake灣中York河口和James河口的水流輸運,發現河口區的水齡垂向分布存在與鹽度類似的層化現象,徑流是控制水齡的主要因素,河口汛期的水齡比平均流量條件下縮短2個月左右。Gong等[5]利用水齡理論,研究了風對Rappahannock河河口環流的影響,發現風對河口水流輸運時間的作用在很大程度上取決于風應力與浮力之間的對比和河口環流的初始狀態。Wang等[14]利用水齡理論探討了長江河口水流輸運時間的洪枯季變化,并系統探討了長江口北槽深水航道工程對河口水流輸運時間的影響。本文集中討論多年平均徑流條件下,潮汐作用對長江河口水流輸運時間的影響。

2 研究區域與研究方法

長江河口水量豐沛,據大通水文站資料,多年平均年徑流總量9 034×108m3(1950—2005年),最大年徑流量13 600×108m3(1954年),最小徑流量6 760×108m3(1978年)。長江徑流量具有明顯的洪、枯季變化,每年5—10月為洪季,徑流占年徑流量的71.7%,11—4月為枯季,徑流量占年徑流量的28.3%。年內最小徑流量出現在1—3月,最大出現在7、8月[15]。

長江河口是中等強度的潮汐河口,口外為正規半日潮,口內為非正規半日淺海潮,平均周期為12時25分。口門多年平均潮差2.66m,最大潮差達4.62m。潮波向口內傳播過程中(傳播方向305°),由于科氏力作用,北岸潮差比南岸大,如南港北岸潮差比南岸高0.4~0.5m。由于口門較寬,納潮量巨大,枯季小潮進潮量為13×108m3,洪季大潮納潮量為53×108m3,兩者相差近4倍。在口外近于平均潮差上游徑流量接近于年平均值時,河口進潮量達266 300m3/s,約為大通站多年平均徑流量的9倍。

本文采用環境水動力學模型(Environmental Fluid Dynamic Code,EFDC)對長江口水齡進行研究。EFDC是在美國國家環保署資助下弗吉尼亞海洋 研 究 所 (VISM,Virginia Institute of Marine Science at the College of William and Mary)的John Hamrick等[16]根據多個數學模型集成開發的綜合模型,該模型集水動力、泥沙輸運、污染物擴散和水質預報于一體,可用于模擬三維動力、物質輸運(包括溫度、鹽度、非黏性和黏性泥沙的輸運)、生態過程和淡水取水、排水等等。其模擬范圍包括河流、湖泊、水庫、濕地、河口以及自近岸到陸架的海域。模式可綜合考慮風、浪、潮、徑流等多因素的影響[13]。

2.1 模型設置

長江河口和杭州灣毗鄰,水沙交換非常頻繁,將它們作為一個整體進行計算不但能更容易給出邊界,而且還可以消除邊界對研究區域的影響,計算區域如圖1,計算網格如圖2。長江口上邊界取在長江中下游最后一個徑流控制站,潮區界安徽大通,距離河口約600km。杭州灣上界取到閘口附近,其余3個外海邊界為東邊界到123°E,北邊界到32°11′N,南邊界到29°53′N。模型區域共有39 727個計算網格,網格的分辨率相差很大,其中最大格距為2 219m(在外海),最小格距僅為254m。在垂直方向分為等間距8層。曲線正交網格在走向上兼顧了長江口的4個入海通道,并兼顧了北槽深水航道工程。這樣一套網格在不過多地增加計算量的同時,卻能保證了關鍵區域的計算精度。時間步長為20s,模型計算90d后,模型達到準穩定條件下,重新啟動,對計算結果進行分析。在數值模擬實驗中,溫度設為常數20℃。

圖1 研究區域與選定斷面

圖2 模型計算區域網格

計算區域的初始條件涉及水位、流速和鹽度,由于水位和流速對外界動力響應較快,初值均取為零;鹽度取洪季或枯季的典型分布值。

外海開邊界通過水位變化作為模式的驅動,水位ζ的表達式為:

式中,A為余水位;fi為各分潮的交點因子;(V0+u)i為分潮的天文相角,它們可由地理位置以及具體的年、月、日求得;ωi為分潮的角頻率;Hi和gi為潮波的振幅和地方遲角;M為分潮個數,在本研究中,M=6,包括S2、M2、N2、K1、O1、P1。

外海開邊界鹽度設為平均流量下全潮平均特征值,其中底層設為恒定值34,表層設為25,通過調節恢復時間Trec使得計算值和實際測量值最接近。

式中,cout為開邊界上一次出流時刻的濃度值,cbnd為指定邊界上的濃度值,tout為由出流轉為入流的時間,Trec為恢復時間[17-18],在該研究中,恢復時間為1h。

本研究主要著重討論長江河口潮汐作用對河口水流輸運時間的影響,長江口上游邊界取多年平均流量(29 000m3/s),鹽度取0,杭州灣上邊界取多年平均流量1 280m3/s,鹽度取0。

河口是由河流進口段、河流河口段和口外海濱組成的一個特殊的水域。就長江河口而言,河流近口段在安徽的大通到江蘇的江陰之間,前者是潮區界,后者為潮流界;河流河口段在江陰到河口攔門沙灘頂之間,也就是滯流點附近。徐六涇在地形上是河口的節點所在,自徐六涇以下河口開始分汊,研究長江口的問題,常以徐六涇作為起點,本文也選擇徐六涇作為水流輸運時間的起點,長江河口水流的輸運時間為相對于徐六涇為起點的時間。

2.2 模型方法

在數學模型中,利用輸運方程計算示蹤物濃度和年齡濃度,考慮示蹤物僅從一個河流邊界進入,不考慮其他源和匯項,因此,示蹤物濃度和年齡濃度可以通過如下方程表示:

式中,C(t,x→)表示蹤物濃度,α(t,x→)為年齡濃度,u為流速,K為擴散系數張量,t表示時間,x→表示空間位置。因此平均年齡可以表示為:

在邊界的設置上,上游邊界示蹤物濃度設為1,年齡濃度設為0,其他的開邊界示蹤物濃度和年齡濃度均設為0。由各層的水齡平均得到垂向平均水齡。

長江河口水動力模型經過大量野外實測數據驗證,具體驗證情況可參考 Wang等[14],結果表明該模型能夠全面反映本海區的主要動力特征,為數值模擬試驗精度提供了保障。選定從徐六涇至河口(122.5°E)通過南港-南槽和南港-北槽的南北兩個縱斷面分析長江河口水流輸運時間的沿程變化(見圖1)。

3 結果與討論

3.1 水齡平面分布

水齡通常用來衡量水體或溶解性物質從進入起到河口某處所需要的時間。圖3給出了多年平均流量下長江河口大小潮全潮平均的垂向平均水齡分布,結果表明,在多年平均流量條件下,水流從進入河口(徐六涇)到南、北港分流口大約需要4d,到最大渾濁帶海域(約122°~122.5°E)需要16d,出長江口(122.5°E)大約需要24d。水齡等值線分布在口門處開始右偏,長江河口水團從徐六涇輸送至杭州灣水域一般需要48d。Wang等[14]研究表明長江下瀉徑流在洪、枯季分別需要20和36d從徐六涇傳輸出長江河口。

長江河口水流輸運在空間分布上具有明顯的橫向結構:受河床摩擦的影響,在大小潮期間均為邊灘水體輸運速度要比主泓處慢。垂向平均水齡在橫向表現為左岸大于右岸(向海方向),水齡等值線右偏,說明水流輸運右岸比左岸更快 (圖3)。同時長江河口水齡空間分布存在與余流分布類似的舌狀結構,水齡則表現為主槽小,相鄰淺灘值大。

圖3 多年平均徑流條件下,長江河口垂向平均水齡的分布

長江河口水齡等值線由徐六涇向下游至口門基本平行分布,水平梯度逐漸變大,在河口口門存在一個水平梯度最大的區域(122°~122.5°E之間),說明此河段水流輸運速度比上、下河段更慢,與長江口最大渾濁帶主體部分吻合,水體輸運速度與最大渾濁帶的形成應該有密切的關系。南匯邊灘位于長江口和杭州灣的分界,潮汐、徑流和灘槽相間的地形決定了本區域復雜的余流結構,水流主槽凈向海,邊灘凈向陸。南匯邊灘處水齡值明顯高于周圍的區域,表明該區域動力條件并不利于水體向外輸運。此外模型結果還表明該區域存在明顯邊灘向陸,主槽向海的順時針環流,前人基于野外觀測數據的分析結果也曾發現類似的結構[19-20]。

3.2 水齡垂向分布

通過分析南北兩個斷面的縱向水齡分布(圖4),發現徐六涇至南北槽分流口附近,水齡垂向分布均勻,水流方向均指向下游;南北槽分流口往下,水齡出現和鹽度相似的分層現象,水齡表層小底層大;水體余流方向表層向海,底層向陸(圖4、圖5)。大潮南北斷面的水齡分層不明顯,垂向上水流方向基本均向海;小潮水齡分層明顯,水流表層向海,底層向陸;北槽上段是底層余流方向轉化的區域,上游部分底層余流向海,下游部分底層余流向陸,即此處為滯留點區域所在,也是攔門沙灘頂。水齡的水平梯度最大的位置也位于此處,說明了本區域水流傳輸速度比上、下河段都要慢。

水流輸運時間的層化與鹽度分層類似,小潮分層更加明顯(圖4)。南北兩個縱斷面水齡分布表明北斷面大潮表底層水齡相差約為2d,其中表層為20d,底層為22d;在小潮期間,表層為16d,而底層為24d,表底層的水齡差異可達到8d。南斷面在大潮期間表底層水齡差異不足2d,而在小潮時段,表底層水齡差異最大可至6d,其中表層水齡18d,底層水齡24d。上述結果表明小潮潮汐混合作用較弱,河口重力環流增強,表層水流向海輸運更快,底層輸運速度則慢。

圖4 多年平均徑流條件下,長江河口水齡的沿程分布(單位:d)

圖5 多年平均徑流條件下鹽度的沿程分布

3.3 垂向平均分布

分析南北縱斷面垂向平均水齡的沿程分布(見圖6),從徐六涇至南北港分流口,水體主要受到淡水控制,且水面寬度變化明顯,在這個河口段水齡與水體輸移的路徑線性增加,對于徑流控制的區域,基本上可用L/ˉu粗略估算水流的輸運時間,其中L為該點與徐六涇的距離,ˉu為斷面平均的余流速。南北港分流口向下游至南北槽分流口區段,由于長興島和橫沙島的影響,水面寬度略微縮窄,水流輸運速度相應增加,該段水齡的梯度略微減小。進入南北槽后,在水面變寬和潮汐混合作用多重影響下,水齡顯著增加,垂向平均水齡的水平梯度變大,在該段水流輸運速度減小,在南北槽下段,水齡增加的速度再次放慢。

圖6 多年平均徑流條件下,長江河口垂向平均水齡的沿程分布

長江河口垂向平均水齡在南北斷面的大小潮過程中表現出相同的變化趨勢。具體表現為多年平均流量條件下,徐六涇至南、北港分流口之間,兩個縱斷面大小潮全潮平均水齡幾乎相等,從南北港分流口向下,大小潮差異越來越明顯,南北斷面,最東端大潮輸移時間比小潮減小3d(圖6)。因此可以認為在南北港分流口以上,長江河口水流輸運時間基本上由河流作用控制,向下開始逐漸受潮汐作用的影響。

從南北兩個縱斷面水齡沿程增加的速度可以發現在距離徐六涇約100km至120km之間,水齡的水平梯度最大,說明該段水流輸運速度比上、下河段都慢,其出現的位置與長江口最大渾濁帶及攔門沙系的主體位置一致。可以認為此段為河口水流輸運時間的徑流主控向潮汐控制的過渡地帶,愈往下游,潮汐作用對水流輸運影響愈加明顯。

3.4 無潮汐作用的水齡垂向分布

為了研究長江河口潮汐作用對水流輸運時間的影響,保持其他條件一致,通過改變開邊界動力邊界條件設置。對比數值計算和實驗結果分析兩種設置下長江河口水流輸運時間變化,確定潮汐作用對該河口水體輸運的作用。不考慮潮汐的振蕩作用時,長江河口水流輸運時間表現出高度的分層,在多年平均流量條件下,表底層水流輸運時間相差高達30d(圖7)。在整個河口段主槽均可以發現表層凈向海,底層凈向陸的垂向環流。距離徐六涇以下80km以下水齡梯度開始顯著增加,水齡等值線非常密集。上述結果表明長江河口的潮汐作用是影響河口水流輸運時間的關鍵要素,河口巨大的進潮量顯著增強河口水流交換能力,影響輸運格局。潮汐混合作用促進上游水體向外海輸運。此外其他的研究也發現長江河口潮致余流非常明顯。Wu等[21]研究表明潮汐的非線性作用與河口淺灘的地形相互影響,潮汐作用加強,潮致余流增加,所以大潮時潮致余流大于小潮。這個結果與斜壓作用占主導作用的弱潮河口不同[4]。

圖7 多年平均徑流,無潮汐作用下長江河口水齡的沿程分布(單位:d)

3.5 無潮汐作用下垂向平均水齡的沿程分布

無潮汐作用下,垂向平均水流輸運時間由潮汐作用下的20d增加至50d,增幅150%。同時,垂向平均水齡的沿程分布也發生了明顯的變化,徐六涇向下20km到80km段水流輸運時間緩慢增加。再向下段,水流輸運顯著減慢,水齡增加速度明顯加快(圖8)。縱斷面自上而下水齡梯度最大的過渡區域在無潮條件并未發現,說明潮汐振蕩作用對最大渾濁帶的形成有重要的影響。通過對比有無潮汐作用的結果發現,當潮汐作用存在的時候,由于潮汐的非線性作用和河口地形相互影響,河口渾濁帶區域水齡梯度大于上下游河段,表明在這個河段,水流輸運速度顯著慢于上下河段,有利于泥沙在此河段堆積,為河口最大渾濁帶形成提供動力支持。上述結果表明:長江河口潮汐的周期性混合作用在調節河口水流輸運時間及物質輸運過程起關鍵的作用。

圖8 多年平均徑流,無潮汐作用下長江河口垂向平均水齡的沿程分布

4 結語

利用長江河口水動力數學模型,研究了多年平均流量條件潮汐作用對長江河口水流輸運時間的影響。研究得出如下主要認識:

(1)在多年平均徑流流量和潮汐作用共同作用下,水流從進入河口到南、北港分流口大約需要4d,到最大渾濁帶海域需要16d,出長江口大約需要24d。南、北港分流口以上,長江河口水流輸運時間基本由河流控制,往下游至最大渾濁帶及攔門沙系存在一個水流輸運速度比上、下游都慢的過渡帶,再向下游潮汐影響比較明顯。

(2)徑流、潮汐共同作用下,長江河口水流輸運時間存在明顯的分層結構,小潮分層比大潮明顯,表底層水齡差異,最大可達6d。

(3)長江河口潮汐的周期性混合作用在調節河口水流輸運時間及物質輸運過程起到了關鍵的作用,無潮條件下水齡表底差異最大可達30d,徑流、潮汐共同作用下水齡水平梯度最大的區域消失,說明潮汐作用與河口最大渾濁帶的形成有密切的關系。

本文利用水流輸運概念量化河口動力過程,成功應用于長江河口。通過數值試驗研究了潮汐作用,可以量化動力條件對河口關鍵過程影響,有助于豐富河口海岸的基本理論。本文研究僅限于多年平均流量與潮汐作用下長江河口水流輸運時間分布特征,風和科氏力作用在大型分汊河口的水流輸運時間作用也有十分重要的意義。

[1] Takeoka H.Fundamental concepts of exchange and transport time scales sea[J].Continental Shelf Research,1984,3:322-326.

[2] Takeoka H.Exchange and transport time scales in the Seto Inland Seas[J].Continental Shelf Research,1984,3:327-341.

[3] Deleersnijder E,Campin J M,Delhez E J M.The concept of age in marine modeling:Ⅰ.Theory and preliminary model results[J].Journal of Marine Systems,2001,28:229-267.

[4] Shen J,Lin J.Modeling study of the influences of tide and stratification on age of water in the tidal James River[J].Estuarine,Coastal and Shelf Science,2006,68(12):101-112.

[5] Gong W,Shen J,Hong B.The influence of wind on the water age in the tidal Rappahannock River[J].Marine Environmental Research,2009,68:203-216.

[6] 匡國瑞,楊殿榮,喻祖祥,等 .海灣水交換研究——乳山東灣環境容量初步探討[J].海洋環境科學,1987,6:13-23.

[7] 朱小兵,高抒,陳妙紅,等 .海南島博鰲港水體交換的初步研究[J].熱帶海洋學報,2003,22:71-77.

[8] 姬厚德,潘偉然,張國榮,等 .筼筜湖納潮量與海水交換時間的計算[J].廈門大學學報:自然科學版,2006,45:660-663.

[9] 董禮先,蘇紀蘭 .象山港水交換數值研究:Ⅱ.模型應用和水交換研究[J].海洋與湖沼,1999,30(5):465-470.

[10] Liu Z,Wei H,Liu G,et al.Simulation of water exchange in Jiaozhou Bay by average residence time approach[J].Estuarine,Coastal and Shelf Science,2004,61:25-35.

[11] 王聰,林軍,陳丕茂,等 .年平均風場作用下大亞灣水交換的數值模擬[J].上海海洋大學學報,2009,18(3):351-358.

[12] Bolin B,Stommel H.On the abyssal circulation of the world ocean:Ⅳ.Origin and rate of circulation of deep ocean water as determined with aid of tracers[J].Deep-Sea Research,1961,8:95-110.

[13] Shen J,Haas L.Calculating age and residence time in the tidal York River using three-dimensional model experiments[J].Estuarine,Coastal and Shelf Science,2004,61:449-461.

[14] Wang Y,Shen J,He Q.A numerical model study of the transport timescale and change of estuarine circulation due to waterway constructions in the Changjiang Estuary,China[J].Journal of Marine Systems,2010,82:154-170.

[15] 陳吉余,沈煥庭,惲才興,等 .緒論[M]//長江河口動力過程和地貌演變 .上海:上海科學技術出版社,1988:1-4.

[16] Hamrick J M.A Three-Dimensional Environmental Fluid Dynamics Computer Code:Theoretical and Computational Aspects[R].Special Report in Applied Marine Science and Ocean Engineering.No.317College of William and Mary,VIMS,1992:63.

[17] Leendertse J J,Gritton E C.A water quality simulation model for well mixed estuaries and coastal seas:Vol.Ⅱ,Computation Procedures[R].Report R-708-NYC,The Rand Corporation(Santa Monica),1971.

[18] Thatcher M L,Harleman R F.A mathematical model for the prediction of unsteady salinity intrusion in estuaries[R].Massachusetts Institute of Technology,Report no.144,1972.

[19] 惲才興 .長江河口潮灘沖淤和灘槽泥沙交換[J].泥沙研究,1983,4:43-52.

[20] 李九發,茅志昌,孫介民 .長江口南匯邊灘及鄰近水域洪季水文泥沙條件分析[M]//長江河口動力過程和地貌演變 .上海:上海海洋學科技出版社,1998:253-267.

猜你喜歡
研究
FMS與YBT相關性的實證研究
2020年國內翻譯研究述評
遼代千人邑研究述論
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
關于遼朝“一國兩制”研究的回顧與思考
EMA伺服控制系統研究
基于聲、光、磁、觸摸多功能控制的研究
電子制作(2018年11期)2018-08-04 03:26:04
新版C-NCAP側面碰撞假人損傷研究
關于反傾銷會計研究的思考
焊接膜層脫落的攻關研究
電子制作(2017年23期)2017-02-02 07:17:19
主站蜘蛛池模板: 欧美日韩精品一区二区视频| 国产成人91精品免费网址在线| 亚洲成人黄色在线观看| 国产精品部在线观看| 欧美色丁香| 婷婷伊人五月| 国产成人AV大片大片在线播放 | 99视频在线观看免费| 色综合网址| 欧美成人h精品网站| 国产不卡网| 99热国产在线精品99| 无码国产偷倩在线播放老年人| 亚洲欧洲日产国产无码AV| 99久久国产自偷自偷免费一区| 久久午夜夜伦鲁鲁片无码免费| 蜜桃视频一区二区三区| 色综合天天操| 成人永久免费A∨一级在线播放| 国产精品亚欧美一区二区三区| 精品91视频| 91外围女在线观看| 久久精品日日躁夜夜躁欧美| 精品人妻系列无码专区久久| 国产精品午夜福利麻豆| 日本精品αv中文字幕| 久草热视频在线| 97成人在线视频| 久久国产V一级毛多内射| 精品欧美视频| 精品色综合| 浮力影院国产第一页| 国产91高跟丝袜| 一级在线毛片| 久久久久国色AV免费观看性色| 国产高清自拍视频| 91免费国产高清观看| 女同久久精品国产99国| 欧美日韩综合网| 四虎影院国产| 亚洲乱码在线视频| 亚洲人免费视频| 亚洲无线国产观看| 91色国产在线| 亚洲二区视频| 91年精品国产福利线观看久久| 欧美一级在线| 2019年国产精品自拍不卡| www.youjizz.com久久| 色亚洲成人| 国产av剧情无码精品色午夜| 草草线在成年免费视频2| 免费看av在线网站网址| 精品视频一区二区三区在线播| 国产高清不卡| 小说区 亚洲 自拍 另类| 制服丝袜一区| 亚洲一道AV无码午夜福利| 国产又粗又爽视频| 福利姬国产精品一区在线| 青青草91视频| 99视频在线免费看| 亚洲国模精品一区| 全裸无码专区| 亚洲男人的天堂网| 91精品视频网站| AV在线天堂进入| 日韩毛片在线播放| 婷婷六月综合网| а∨天堂一区中文字幕| 综合网天天| 亚洲中文字幕手机在线第一页| 色综合天天娱乐综合网| 91精品国产一区自在线拍| 无码国产伊人| 亚洲成aⅴ人片在线影院八| 九九九久久国产精品| 极品尤物av美乳在线观看| 欧美综合中文字幕久久| 九一九色国产| 中文字幕第4页| 97超碰精品成人国产|