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

高黎貢山復雜地形下局地環流的數值模擬研究

2021-08-06 02:24:14邢雯雯孫績華劉輝志許魯君
大氣科學 2021年4期

邢雯雯 孫績華 劉輝志 許魯君

1 廣東省氣象臺,廣州 510080

2 中國科學院大氣物理研究所大氣邊界層物理和大氣化學國家重點實驗室,北京 100029

3 云南省氣象臺,昆明 650034

4 中國科學院大學,北京 100049

1 引言

局地環流的形成與演變不僅受到熱力和動力作用的影響,還與下墊面的植被、環境風等有關(Arritt and Pielke,1986;Gross,1987;Barr and Orgill,1989;Clements el al.,1989;Doran,1991)。局地環流系統對有組織的對流系統的生成、局地氣候、空氣污染的傳輸與擴散、霧的形成與消散、低云的形成等有重要作用。研究局地環流和相關大氣邊界層的特性對于氣象預報和空氣質量的評估是很重要的(Prasad et al.,2017)。

山脈對氣流不僅有顯而易見的阻擋作用,還會在背風坡產生渦旋氣流、重力內波(劉輝志等,2001)、波破碎(齊瑛和傅抱璞,1992)等動力強迫現象,還會形成山谷風、坡風、下泄流、梵風等等。早在上世紀50年代,葉篤正先生(1956)在小地形對于氣流的影響的綜合報告中總結到山地流場因風速不同會形成四種結構:無擾動的氣流、渦旋氣流、波狀氣流和轉子氣流。比起平原的山谷,高原山谷的山谷風環流更加強盛(Schmidli and Rotunno,2012);還有封閉型山谷比開放型山谷的山谷風環流發展厚度大且特征顯著(王瑾等,2012)。這些研究多利用數值模擬或觀測分析揭示地形對氣流的影響。然而,局地風的發展演變還會受到環境風場的影響。眾所周知,在弱背景天氣條件下,更有利于局地環流的發展。在阿爾卑斯山脈的因河河谷,在偏南風背景下山谷風會抵抗環境風,且氣流在阿爾卑斯山南面上升,北面下沉(Z?ngl,2009)。由背景風和山谷風的風向關系還可以了解山谷風的驅動機制(Cogliati and Mazzeo,2006)。

WRF(weather research and forecasting)模式被廣泛運用于研究山谷地形中大氣邊界層的風和湍流特征以及大氣邊界層垂直結構。在這些研究中,由于地域、氣候、地形等差異,效果較好的參數化方案也各不相同(劉振鑫等,2012; 屠妮妮等,2012;張小培和銀燕, 2013;尹相玉,2014;李斐等,2017)。我國地形西高東低,橫斷山脈位于青藏高原東南緣,平均海拔4000米以上,不僅地表粗糙度大,在動力上對大尺度環流有阻擋作用,而且熱力上隨季節變化和日變化差異都很大,這就導致了山脈東西部的氣候差異以及山區氣候明顯的垂直變化。高黎貢山這種位于西南季風和東南季風交匯區的南北走向的狹長型山脈,其附近的局地環流特征研究目前開展得不多。

高黎貢山是橫斷山脈最西部的山脈,北至西藏自治區,南達云南保山市境內,南北長約135千米,平均海拔3500米,南段海拔稍低,東西僅9千米。其西側的深切型怒江大峽谷最低海拔700米,相對高差大。該區域位于亞洲季風前沿,地處東亞水汽通道上游,受高空西風環流、印度洋和太平洋季風環流的影響。孟加拉灣西南季風對該區的影響占主導地位(李汀和琚建華,2013)。已有研究指出,濕季除西南季風外,東南季風也會影響云南部分地區(陳艷等,2006)。而影響云南的太平洋東南季風也基本分為三支:西太平洋副熱帶高壓南側的東南氣流、南海赤道輻合帶(ITCZ)北側的東南風以及臺風外圍氣流。但只有西太平洋副熱帶高壓伸過105°E時,東南季風和攜帶的水汽才由副高南部外圍從西太平洋和南海進入云南(徐嘉行等,1984)。騰沖站和保山站被分割于高黎貢山東西兩側,受背景風和地形的影響,保山與騰沖兩地的氣溫、相對濕度、降水、局地環流都存在差異,形成了不同的局地氣候。

本文利用騰沖氣象站的地面觀測資料、探空資料和保山氣象站的地面觀測資料,分析兩站的氣候條件差異;結合中尺度數值模式,對比不同邊界層參數化方案(YSU、MYJ、MYNN3、ACM2、Boulac)的模擬效果,研究WRF模式在騰沖保山地區復雜地形下的適用性,并選出最佳方案對干季和濕季的局地環流進行數值模擬,分析不同環流背景下的局地環流特征。

2 數據來源、試驗設計及邊界層參數化方案簡介

2.1 數據來源

本文采用高黎貢山兩側的騰沖國家基準氣候站(下稱騰沖站)及保山市氣象站(下稱保山站)2016年的資料進行分析及模式結果驗證。騰沖站(24.983°N,98.5°E)位于云南保山市。騰沖站海拔1697 m,東距高黎貢山山脈約20 km,場地四周為騰沖市郊區農田,所測資料基本代表了高黎貢山以西的邊界層特征。保山站(25.117°N,99.183°E)位于高黎貢山以東,海拔1652.2 m,西距高黎貢山約32 km。兩個站主要是常規氣象要素觀測。騰沖站與保山站的位置如圖1b所示。

2.2 試驗設計

本文利用WRF3.7(Weather Research and Forecasting Model v3.7)模式,對比分析YSU、MYJ、MYNN3、ACM2、BouLac五種邊界層參數化方案在高黎貢山的適用性,模擬時段為2016年12月3日08時(北京時,下同)至8日08時;并采用最優邊界層參數化方案對2016年高黎貢山干季和濕季局地環流特征進行模擬分析。干季模擬時段同上,濕季模擬時段為2016年9月23日08時至28日08時。模式采用四層嵌套,最內層100 km×100 km,分辨率1 km,中心點位于高黎貢山南段山頂。垂直分層35層,2 km以下10層。模式輸入的初始場采用ECMWF ERA-interim的0.75°×0.75°再分析資料,地形資料為Modis資料,并且采用了次網格地形參數化方案(Jiménez and Dudhia,2012)來提高模擬區域內的近地面風的模擬 效 果(Jiménez and Dudhia,2013;Lee et al.,2015)。模式的四層模擬中采用的物理方案是一致的,短波輻射方案均采用Dudhia方案(Dudhia,1989),長波輻射方案采用RRTM方案(Mlawer et al., 1997),微物理方案采用WDM6方案(Lim and Hong,2010),積云對流方案為Kain-Fritsch方案(Kain and Fritsch,1993),陸面過程采用Noah方案(Chen and Dudhia,2001)。圖1a為模擬區域。圖1b為最內層模擬區域的地形高度,包括了高黎貢山的南段。中間海拔3200 m以上的山脈是高黎貢山,其東側怒江大峽谷最低海拔低于800 m,峽谷東側是碧羅雪山。地勢呈現北高南低,山谷相間分布。圖1c是地表利用類型分布圖。在中西部的山脈地區主要是闊葉林、針葉林和混合森林,其中也混有一部分的草地和灌木。東部主要是草地和灌木,城市周邊分布著農田,西北地區也有農田分布。植被會影響地表能量通量的變化以及大氣邊界層內的其他變量。

圖1 (a)四重模擬區域;(b)d04區域模擬地形高度(陰影)及站點位置(GLGM代表高黎貢山,T代表騰沖站,B代表保山站);(c)d04區域MODIS20類土地利用情況Fig.1(a)Model domains(d01,d02,d03,and d04);(b)terrain height of d04 and locations of two stations(GLGM is Gaoligong mountains,T is Tengchong station and B is Baoshan station);and (c)land usetype of d04 from MODIS 20-classdataset

分析時段內的天氣形勢如圖2所示,12月4日850 hPa等高線稀疏,云南地區沒有明顯的影響系統。海平面氣壓場顯示四川,貴州東部地區存在地面高壓。從9月27日的環流形勢來看,850 hPa我國北方地區受高壓控制,華南沿海有熱帶氣旋正在登陸,低壓中心將逐漸向西移動。

圖2 2016年12月4日08時(北京時,下同)(a)850 hPa位勢高度場(單位:dagpm)和(b)海平面氣壓場(單位:hPa);2016年9月27日08時850 hPa(c)位勢高度場(單位:dagpm)和(d)海平面氣壓場(單位:hPa)(G代表高中心,D代表低中心)Fig.2(a)850-hPa geopotential height field(units:dagpm)and(b)sea-level pressure field(units:hPa)at 0800 BJT 4 December 2016,(c)850-hPa geopotential height field (units:dagpm)and (d)sea-level pressurefield (units:hPa)at 0800 BJT 27 September 2016.G standsfor high center,D stands for low center

2.3 邊界層參數化方案簡介

不同的邊界層參數化方案會導致邊界層內的溫度、濕度以及降水模擬結果存在明顯差異(Xie et al.,2012)。本文選用常用的YSU、MYJ、MYNN3、ACM2、BouLac五種方案進行對比,分析模擬效果,選出最佳方案。

YSU方案(Hong et al.,2006)是基于K理論的一階非局地閉合方案。YSU方案解決了MRF方案中大氣層結過于穩定的問題。該方案顯式處理夾卷過程,降低了動力強迫的湍流混合而增加了熱力強迫誘導的湍流混合。ACM2方案(Pleim,2007)是一階閉合非局地方案,在ACM1方案的基礎上增加了局地湍流輸送,可以模擬由浮力作用引起的非局地向上輸送和局地湍流交換。MYJ、MYNN3、BouLac方案是基于湍流動能TKE的方案。有研究指出,TKE方案比一階閉合方案更能反映邊界層的整體湍流結構(Holt and Raman,1988)。其中MYJ方案(Mellor and Yamada,1982)通過迭代求解TKE產生和損耗微分方差計算邊界層高度。MYNN3方案(Nakanishi and Niino,2004)加入了凝結物理過程,改進了湍流動能和混合層高度偏低的問題。BouLac方案(Bougeault and Lacarrère,1989)根據湍流動能的垂直分布計算邊界層高度。

3 結果分析

3.1 高黎貢山東、西兩側氣象特征

對比分析2016年高黎貢山東、西兩側騰沖和保山站地面氣象觀測資料,發現全年氣溫8月最高,1月最低;最高溫度出現在14時,最低氣溫出現在7時;兩站的溫度年變化和日變化趨勢相似,但保山站的日平均氣溫明顯高于騰沖站,平均高出1.9°C(表1),年平均的溫度日變化保山也明顯高于騰沖(圖3b)。兩站的相對濕度全年50%~100%,其中3月相對濕度最??;并且保山站的相對濕度低于騰沖站,平均低了約14%(表1),相對濕度日變化保山也低于騰沖(圖4b)。

圖3 2016年騰沖站(黑色)和保山站(灰色)氣溫(單位:℃)(a)年變化和(b)日變化Fig.3(a)Annual variation and (b)daily variation of air temperature (unit:°C)at Tengchong station (black)and Baoshan station (gray)in 2016

圖4 2016年騰沖站(黑色)和保山站(灰色)相對濕度(a)年變化和(b)日變化Fig.4 (a)Annual variation and (b) daily variation of relative humidity at Tengchong station (black)and Baoshan station (gray) in 2016

表1 騰沖站和保山站全年、干季和濕季的日平均氣溫(單位:°C)和相對濕度Table 1 Daily average air temperature(units:°C)and relative humidity at Tengchong and Baoshan stations throughout the year,dry and wet seasons

騰沖站全年平均風速1.94 m s?1,干季(圖5a)風向以東南風為主,濕季有西南風也有東南風。風向日變化表現出上午有短時的西北風,白天多為西南風,夜間轉為東南風;保山站全年平均風速1.75 m s?1,干季(圖5c)主要吹西南風,濕季為東南偏南風、偏南風和東南風。無論騰沖還是保山,西南風風速多為4~7 m s?1,大于東南風主要風速2~4 m s?1。

圖5 2016年騰沖站(a)干季、(b)濕季和保山站(c)干季、(d)濕季風玫瑰圖Fig.5 Wind rosesof (a) dry season and (b)wet season at Tengchong station and (c)dry season and (d) wet season at Baoshan station in 2016

南北走向的高黎貢山阻礙了西南季風,其中一部分氣流被阻擋,一部分較強氣流翻越過高黎貢山。來自孟加拉灣的西南季風,在迎風坡爬升時水汽凝結降落,氣流過山以后溫度升高而濕度降低。所以位于山脈西側的騰沖站相比于東側的保山站溫度低且濕度高。哀牢山同樣在濕季西坡濕度資源明顯優于東坡(劉玉洪等,1996)。

3.2 不同邊界層參數化方案的適用性

為了得到最佳的模擬效果,本文選擇的5種行星邊界層方案對比的結果如下。溫度模擬的誤差主要是最低溫度模擬偏高,YSU方案和ACM2方案模擬的平均絕對誤差最小(表2)。觀測風向受地理位置和局地環流影響很大。上午10~11時常出現西北風,白天多為西南風,夜間轉為東南風。模擬結果顯示,MYNN3模擬的風向平均絕對誤差最小,MYJ和BouLac方案其次,但是YSU和MYJ方案的日變化趨勢與觀測值更加接近。騰沖站觀測平均風速1.8 m s?1,各個方案模擬的平均絕對誤差都偏小,其中ACM2最小。綜合風速樣本結構分析,MYJ方案和ACM2方案模擬效果較好(圖6d)。

圖6 2016年12月4日觀測值與5種邊界層參數化方案模擬值對比:(a)2 m溫度(單位:°C);(b)10 m風向[單位:(°)];(c)10 m風速日變化(單位:m s?1);(d)10 m風速箱圖Fig.6 Comparison of observed values on 4 December 2016 with the simulated values of five boundary layer parameterization schemes:(a)2-m temperature (units:°C);(b)10-m wind direction [units:(°)]of daily variation;(c)10-m wind speed (units: m s?1);(d) box chart of 10-m wind speed

表2 五種邊界層參數化方案對溫度、風向和風速的模擬值與觀測值的平均絕對誤差Table 2 Mean absolute errors between observed and simulated values of air temperature,wind direction,and wind speed by the five boundary layer parameterization schemes

3.3 干季和濕季局地環流分析

3.3.1干季局地環流

基于3.2節的分析,干季和濕季局地環流的模擬采用了YSU邊界層參數化方案。最內層模擬區域內的2 m溫度疊加10 m風矢量圖(如圖7),黑色等值線為地形高度。夜間(圖7a)山體溫度較低,高黎貢山山頂氣溫低至零下,其東側的怒江大峽谷溫度較高,山谷之間溫差較大。騰沖和保山(圖7d中黑色圓點所示位置)由于城市熱島效應溫度偏高。夜間風速較小,山坡輻射冷卻降溫較大,山谷降溫較慢,盛行由山坡吹向山谷的山風,高黎貢山東西兩側山風明顯,西北部海拔較高的區域和東部的碧羅雪山,也以山坡吹向山谷的山風為主要環流形勢,山頂氣流輻散而山谷氣流輻合。并且,東側怒江大峽谷內風向偏北,直至上午模擬區域內的氣流基本都有偏北分量,日出之后地面開始接收太陽輻射加熱,山谷開始升溫,09時山風減弱而谷風開始發展。山坡東側升溫快,在氣壓梯度力的作用下,高黎貢山東側最先開始形成由山谷吹向山坡的谷風。隨著溫度升高,10時(圖7b)西側的谷風也逐漸發展起來。14時(圖7c)谷風環流更加明顯,高黎貢山兩側氣流沿山坡上升,氣流在山頂附近輻合,怒江大峽谷氣流輻散。騰沖站基本為西南氣流,保山站為偏東氣流。高黎貢山西側氣流轉為有偏南分量。19時高黎貢山東坡開始出現山風環流,20時西坡也開始出現山風,21時山風環流加強。夜間高黎貢山山體氣流輻散,騰沖站主要為南風,東側保山站主要為西北風。綜上,地面風速較小,尤其在夜間,高黎貢山東西兩側主要在地形熱力作用下形成有日變化的山谷風環流;從夜間到白天,氣流的偏北分量逐漸轉為偏南分量,同樣是北高南低的地形熱力作用導致。夜間山頂氣流輻散,山谷氣流輻合,白天則相反。騰沖站和保山站分別位于山的西側和東側谷內,風向受山谷風環流發展轉換的影響,風向基本相反,這與觀測事實是相符的。

圖7 模擬的2016年干季(12月4日)10 m風場(箭頭,單位:m s?1)、2 m溫度(陰影,單位:°C)和地形高度(等值線,單位:m):(a)02:00;(b)10:00;(c)14:00;(d)20:00Fig.7 Simulated the 10-m wind field(black vector,units:m s?1),2-m temperature(shaded,units:°C),and terrain height (contour,units: m)in 2016 dry season (December 4th):(a)0200 BJT;(b)1000 BJT;(c) 1400 BJT;(d)2000 BJT

為了分析位溫和垂直環流結構,沿25.23°N做5 km以下的垂直剖面圖(圖8見文后彩圖)。干季中層盛行偏西風。夜間,低層垂直剖面圖上全部為平直的西風氣流。在高黎貢山東側的怒江大峽谷,西風過山以后形成了小的空腔區和逆向氣流。西側山谷低層有被高黎貢山阻擋的折返氣流。日出以后,對流邊界層開始發展,在高黎貢山西側山谷地面以上1 km出現東風氣流,西側怒江大峽谷內的逆向氣流也更加旺盛。由于東側碧羅雪山的阻擋,峽谷內近地層位溫較低,上層背風坡升溫快,高黎貢山兩側等位溫線東高西低。谷風環流也開始發展。14時,西風更加強盛,且在西側山谷內氣流上升往兩側下沉,形成小的渦旋,對流邊界層高度達1 km以上。高黎貢山東側大峽谷內西風氣流逆轉下沉,碧羅雪山西坡出現下坡風。日落后地面長波輻射冷卻,谷風減弱而山風逐漸發展。4日晚上冷高壓南壓,東北季風較強,20時碧羅雪山以西3.2 km以下出現東風氣流,東西風交匯,湍流混合層較弱。峽谷內東風下傳,西側山谷內氣流仍然存在小的渦旋。

圖8 2016年干季(12月4日)沿25.23°N的位溫(填色,單位:K)、風場(箭頭,單位:m/s)的垂直剖面(黑色陰影是海拔超過800 m的地形):(a)02:00;(b)10:00;(c)14:00;(d)20:00Fig.8 Vertical cross sections of potential temperature(shaded,unit:K)and wind field (vector, unit: m s?1)at 25.23°N in 2016 dry season (December 4th)(Black shadows are terrain over 800 m above sea level):(a)0200 BJT;(b)1000 BJT;(c)1400 BJT;(d)2000 BJT

3.3.2濕季局地環流

濕季騰沖保山地區受偏東背景風的影響,夜間高黎貢山山頂已轉為偏東風,在東風的影響下東側的山風明顯比西側弱(圖9a見文后彩圖)。與干季不同的是西側仍以山風環流為主而風向偏北,但東側風向偏南。白天,谷風環流10時發展旺盛,谷風在山頂輻合。14時,高黎貢山東側的偏東谷風加大,甚至在有些地方越過山頂,明顯抑制了西側偏西谷風的發展(圖9c見文后彩圖)。在青藏高原的南克湖盆地,較強的背景風會推遲湖風甚至模擬不出湖風(Gerken et al.,2014)。騰沖站風向西南,保山站風向東南。綜上,東風導致白天高黎貢山西側谷風減弱,夜間東側側山風減弱??梢姰敱尘帮L較強和風向與局地風風向相反時,會抑制局地環流的發展。

圖9 同圖7 ,但為2016年濕季(9月27日)的10 m風場、2 m溫度和地形高度Fig.9 Sameas Fig.7, but for 10-m wind field,2-m temperature,and terrain height in 2016 wet season (September 27th)

夜間4.4 km到地面全部為東風控制。上層西風與下層東風交匯形成交換帶(圖10見文后彩圖)。上午10時山谷風環流開始形成。東西風相遇的高度下降至3.5 km,在東、西風相遇高度上形成若干個小的渦旋。在東側峽谷內東風受高黎貢山阻擋后下沉或彎折。午后高黎貢山西側高層西風動量下傳,與下層東風產生擾動,西側近地面氣流起伏波動,等位溫線也隨之起伏。東側的東西風交換帶更低一些,東風在怒江大峽谷下沉,怒江大峽谷以314 K位溫充分發展混合層。東風氣流在低層被高黎貢山阻擋而高層強西風氣流翻越過高黎貢山逐漸下傳,使得山脈兩側東、西風混合層高度不同。20時氣流已沒有中午時那么強盛,東風高度在4 km左右。濕季,在偏東背景風的影響下,中層東風和西風的動量交換會使底層產生波動,東風氣流的影響高度僅在5 km以下。

3.4 邊界層高度分析

圖11為干季和濕季騰沖站、保山站邊界層高度的日變化。夜間穩定邊界層高度幾乎保持不變,維持在幾十米的高度。09時谷風環流開始時邊界層開始發展,午后對流邊界層充分發展,14時至15時高度達到最高,干季騰沖站大氣邊界層高度可達到2500 m以上,保山站大氣邊界層高度為2169 m,騰沖邊界層高度比保山高300 m左右。干季邊界層高度高于濕季。騰沖站和保山站都位于高黎貢山南段緯度上,基于上述干季和濕季位溫和風場的剖面分析,可以注意到干季背景風環流較弱,利于局地環流的發展;而濕季受偏東背景風影響時,白天的谷風環流與背景風風向是相反的,因而抑制了局地環流的發展,邊界層高度也就低于干季。正如當背景風風速較大時,湖風會被推遲或者抑制(Gerken et al.,2014)。另外,高黎貢山西側干季由于山脈阻擋,西風氣流下沉會形成小的渦旋,湍流混合充分,邊界層高度高于保山;濕季受偏東風影響時,背景風使高黎貢山西側谷風和東側山風減弱,而湍流混合層主要受白天谷風發展的影響,故而濕季騰沖與保山的邊界層高度相差不大。

圖11 模式中干季(2016年12月4日)(黑色)、濕季(2016年9月27日)(灰色)騰沖站(方框)和保山站(圓點)的邊界層高度(單位:m)對比Fig.11 Comparison of boundary layer height(units: m)between Tengchong station(square)and Baoshan station(dot)of dry season(4 December 2016)(black)and wet season(27 September 2016)(gray)in model

4 結論

本文基于高黎貢山東、西兩側騰沖和保山站氣象觀測資料,驗證WRF模式不同邊界層參數化方案在復雜地形下的適用性,并模擬分析高黎貢山局地環流特征。研究表明,高黎貢山西側比東側溫度低、濕度大,東、西兩側干季和濕季主導風向不同。高黎貢山南段的騰沖—保山地區谷風出現在09時,山風出現在19時。局地環流以山谷風環流為主,白天多為偏南風,夜間為偏北風和偏西風。白天山頂氣流輻合,山谷輻散,夜間相反。白天風速大于夜間。干季西風風力較弱,有利于對流層低層局地環流的發展。濕季偏東風較強,谷風環流發展受到抑制,邊界層高度低于干季。另外,干季較強的西風遇到高黎貢山,在西坡下沉和形成渦旋,西側湍流混合充分,邊界層高度高于東側。濕季受較強的偏東背景風削弱了西側谷風環流的發展,因而騰沖與保山的邊界層高度相差不多。

本文對高黎貢山南段地區的局地環流做了個例分析,結果與該區的氣候特征、背景風等較為一致,具有一定的代表性。但是僅模擬了兩個時段的局地環流,若要了解高原復雜地形下在西南季風和東南季風交匯區的局地環流,還需要匯集更多的地面觀測資料以及背景風資料來進行更加深入的研究。

主站蜘蛛池模板: 国产菊爆视频在线观看| 久久中文字幕2021精品| 国产亚洲精久久久久久久91| 99久久精品国产综合婷婷| 尤物在线观看乱码| 六月婷婷精品视频在线观看| 亚洲专区一区二区在线观看| 亚洲,国产,日韩,综合一区| 国产v精品成人免费视频71pao| 日本黄网在线观看| 欧美日韩亚洲国产主播第一区| 亚洲日本在线免费观看| 日韩AV手机在线观看蜜芽| 欧美在线免费| 色婷婷久久| 青青草久久伊人| www亚洲精品| 亚洲中文字幕无码mv| 婷婷色婷婷| 手机成人午夜在线视频| 综合久久五月天| 国产噜噜噜| 熟女日韩精品2区| 欧美性猛交一区二区三区| 成人午夜视频网站| 欧美在线天堂| jizz在线免费播放| 影音先锋亚洲无码| 91探花在线观看国产最新| 欧美全免费aaaaaa特黄在线| 国产精品无码AⅤ在线观看播放| 丝袜亚洲综合| 激情在线网| 国产区91| 好紧好深好大乳无码中文字幕| 77777亚洲午夜久久多人| 本亚洲精品网站| 露脸真实国语乱在线观看| 五月激情婷婷综合| 精品偷拍一区二区| 影音先锋丝袜制服| 99视频在线精品免费观看6| 国产成人综合在线观看| 国产精品密蕾丝视频| 美女毛片在线| 青青草原国产精品啪啪视频| 激情网址在线观看| 精品久久高清| 国产xxxxx免费视频| 亚洲资源站av无码网址| 亚洲成在人线av品善网好看| 91精品视频网站| 亚洲精品成人7777在线观看| WWW丫丫国产成人精品| 国产精品9| 国产九九精品视频| 伊人查蕉在线观看国产精品| 国产精品无码作爱| 亚洲福利片无码最新在线播放| 亚洲精品第一在线观看视频| 精品一区二区无码av| 欧美午夜一区| 亚洲精品在线91| 免费AV在线播放观看18禁强制| 国产无码高清视频不卡| 一本大道无码高清| 国产成人超碰无码| 免费av一区二区三区在线| 国产区精品高清在线观看| 国产美女一级毛片| 国产精品女主播| 亚洲专区一区二区在线观看| 69精品在线观看| 91网址在线播放| 亚洲va在线∨a天堂va欧美va| 精品91视频| 天堂岛国av无码免费无禁网站 | 国产日韩精品一区在线不卡| 国产黑丝一区| 亚洲天堂2014| 夜夜爽免费视频| 亚洲美女一区|