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

垂直湍流輸送對大洋的重力位能和混合過程的影響

2014-08-04 00:51:16張宇1林一驊2王輝1
大氣科學 2014年5期

張宇1 林一驊2 王輝1

?

垂直湍流輸送對大洋的重力位能和混合過程的影響

張宇林一驊王輝

1國家海洋環境預報中心,北京100081;2中國科學院大氣物理研究所大氣科學和地球流體力學數值模擬國家重點實驗室,北京100029

利用WOA09(World Ocean Atlas 2009)全球大洋溫鹽客觀分析數據,計算了不同湍流垂直混合系數下全球大洋重力位能的變化,并分析了混合系數、浮力頻率和重力位能變化之間的關系。在此基礎上,進一步探討了湍流混合造成的能量轉換對湍流參數化的影響。結果表明,大洋中的垂直湍流運動不僅僅是動能能匯,而且是一個重要的外部能量轉化為重力位能的途徑。垂直湍流增加的重力位能在混合系數取0.1 cms時為0.08 TW,參考前人研究結果,外部能量輸入甚至可引起等效于全球平均12 cms的垂直混合系數。一般而言,層結越穩定、混合系數越大,垂直湍流對重力位能的影響也越大。考慮湍流動能可轉化為重力位能后,參數化方案可以得到和實際觀測更接近的湍流動能耗散率和混合系數。

混合系數 重力位能 湍流動能耗散率 湍流混合參數化

1 引言

垂直湍流混合是海洋中最重要的物理現象之一,也是大洋環流模式需要解決的關鍵問題。在現有數值模式的分辨率條件下,解決湍流混合問題只能依靠次網格參數化方案。在各種參數化方案中,垂直湍流混合系數(下文簡稱混合系數)的時空分布成為最關鍵也是最難解決的科學問題。長期以來,不同學者從理論和觀測角度嘗試解決這一難題。Munk and Wunsch(1998)和Munk(1966)從大洋中的平均垂直速度和維持溫鹽分布的大尺度平衡需要出發,推算出了全球大洋平均混合系數量級應為1 cms。但是后來的研究尤其是觀測事實表明(Gregg,1989;Ledwell et al., 1998;Osborn,1980),在大洋內部混合系數量級普遍為0.1 cms,比大尺度平衡理論的結果小了一個量級。這說明垂直湍流混合強度在大洋中分布極不均勻,必然存在混合系數超過1 cms的強混合區域。

由于現場直接觀測混合系數技術困難、代價昂貴,一直以來只有極少量的直接觀測數據,遠不能滿足大洋環流的研究需要。因此,許多學者發展了從其他較易觀測到的物理量中推算出混合系數的方法。Osborn(1980)的經典工作使得從湍流動能耗散的觀測數據中推算出混合系數成為可能,極大推動了實際大洋中垂直湍流混合的觀測研究。St Laurent et al.(2012)、Nagasawa et al.(2007)、Thompson et al.(2007)等都采用Osborn(1980)的方法結合觀測數據得出了大洋中特定區域混合系數的時空分布。

傳統意義上認為湍流運動主要起能量耗散的作用,大量的數值模式僅從考慮海水運動的平均狀態誘發湍流難易程度的角度來參數化湍流運動,而不考慮外部能量輸入對湍流的影響。風和潮汐是大洋最主要的外部能量來源,在風能和潮汐能輸入強的地區應當存在較強的垂直湍流,這個觀點也已經被觀測和理論研究所證實(Ledwell et al., 2000; Wang and Huang,2004a,2004b)。Craig and Banner (1994)由此從Mellor and Yamada(1982)的湍流混合理論出發,提出了從海洋表層風和波浪特征推算大洋中風浪引起湍流動能垂直耗散的參數化方案,下文簡稱CB94。CB94和Osborn方案相結合,可以建立由風、海浪估算海洋上層混合系數的一整套參數化方法。但是CB94方案并未考慮到海水密度的變化,因而也無法描述由密度變化而引起的重力位能改變,這在CB94方案所關注的海洋上層100 m垂直范圍內是近似合理的。

事實上,垂直湍流混合不僅會耗散能量、混合溫鹽,還可以通過改變大洋中的溫鹽層結進一步影響大洋中海水的重力位能(GPE:Gravitational Potential Energy)。一般情況下,由于大洋中以穩定層結為主,垂直混合會增加海水重力位能,因此大洋中的垂直混合不僅僅是動能能匯,也是重力位能能源(黃瑞新,1998)。在海洋中風浪較強的區域如南極繞極流(ACC:Antarctic Circumpolar Current),風浪可以影響到1000 m水深處的湍流混合,在這個深度上混合對重力位能的貢獻有可能改變湍流動能耗散強度的垂直分布,從而進一步影響參數化的準確性。因此,本文將從能量平衡的角度出發,細致分析垂直混合對海水重力位能的影響,并進一步研究其對類似CB94的方案參數化100 m至風浪影響最大水深范圍內混合系數的影響。

2 資料與計算方法

2.1 資料來源

本文使用了由美國國家海洋數據中心(NODC:National Oceanographic Data Center)提供的WOA09(World Ocean Atlas 2009)海洋溫鹽客觀分析數據。該資料的覆蓋范圍為(89.5°S~89.5°N,0.5°E~0.5°W),水平分辨率為1°×1°,垂向0~5500 m分為33層,包含現場溫度,實用鹽度,溶解氧等數據,時間分辨率為氣候態年平均和月平均。此外,還使用了由歐洲中期天氣預報中心(ECMWF:The European Centre for Medium-Range Weather Forecasts)提供的ERA Interim(ECMWF Re- analyses Interim data)資料中的海表風應力和海表面波浪高度數據。數據空間分辨率為1.5°×1.5°,時間分辨率為6小時。

為滿足本文研究需要,對原始數據做了預處理工作,將垂直上的不穩定層結通過垂直對流調整為穩定層結。本文用The Gibbs Sea Water(GSW)Oceanographic Toolbox of TEOS-10標準程序將WOA09中的現場溫度轉換成位溫(以海平面為參考面),實用鹽度轉換為絕對鹽度,并計算了每一層的現場密度和以海平面為參考面的位密,用于進一步計算海洋上層的浮力頻率。海洋中深層則用每一層深度為參考面計算上下相鄰兩層的位密,并進一步計算浮力頻率。

圖1 垂直混合對重力位能改變貢獻的計算方法圖示

2.2 計算方法

本文的計算主要為混合系數引起的重力位能改變和風輸入海洋上層的湍流動能耗散率及混合系數。

混合系數引起重力位能的改變計算流程如下:

(1)計算初始狀態海水的重力位能,以WOA09中的海底最深處為參考面(5500 m)。

式中,,,分別對應經向、緯向和垂直方向上的網格點編號,,,對應經向、緯向和垂直方向上的最大格點數。為大洋海水總重力位能,m,,為格點質量,為重力加速度,本文取9.8 m s。H,,為格點質心和5500 m水深間的距離。

(2)計算混合后格點新的位溫、鹽度和現場密度。

式中,,,,,分別為位溫、絕對鹽度、混合系數、壓力和海水現場密度,上標?表示經過混合后的新物理量,表示混合所用時間,本文取為1 s。對真實微尺度湍流而言,這個時間尺度較大。但在大洋環流數值模式中,這個尺度相當于積分步長,由于溫鹽湍流理論參數化過程普遍在模式的斜壓模態中向前積分,而模式的斜壓積分步長相對較長,在水平網格分辨率為1°,垂直分辨率為1~10 m時,斜壓積分步長量級為10s,因此取值1s是合適的。此外,取值1s也便于下文計算結果直接以功率單位(W)表達。

(3)計算混合后格點的新高度?和大洋新的重力位能?及重力位能改變。

(4)除了大洋總重力位能的變化,我們還需要詳細了解每一層垂直混合對重力位能改變所做的貢獻。實際計算中混合總是位于格點的界面上,因此每一界面的混合可細分為將下層格點的物質混合至上層及將上層格點物質混合致下層兩個過程(圖1)。如果將這兩個過程分開考慮,那么經混合后上、下兩層格點新的位溫應為

(4)

式中,下標up和down分別代表混合界面上、下格點的物理量,為格點厚度。類似的,我們也可導出上、下層相應的、、、以及相應的上、下層格點厚度因混合的變化、。每一格點重力位能的貢獻可以視作格點厚度的改變推動格點上方水柱上下位移而做的功。因此,每一層混合對整個大洋的重力位能貢獻為:

風輸入海洋上層的湍流動能耗散率計算方案主要參照Baumert(2005)的方案,將海洋上層分為波浪混合層、湍流混合層和壁面律層(該層波浪混合效應按照普朗特湍流假設處理)分別計算(圖5a)。波浪混合層計算公式為

湍流混合層計算公式為

(7)

3 常混合系數對大洋重力位能的影響

本文首先計算了常混合系數下,全球大洋重力位能的變化,結果見圖2、圖3、圖4。從圖2中可看出,在常混合系數取0.1 cms時,垂直混合引起的全球大洋重力位能分布非常不均勻,熱帶地區最高,北冰洋沿岸次之,南大洋地區最低。孟加拉灣、巴拿馬灣、喀拉海和拉普捷夫海等有河流大量淡水注入的海域重力位能改變最顯著。南中國海、東印度洋、赤道西太平洋等有大量熱量輸入的地區重力位能改變也較為顯著。這種和溫鹽密切相關的空間分布形態說明了常混合系數下垂直湍流混合對重力位能的貢獻受海洋層結的直接影響。全球平均的混合對重力位能貢獻垂直分布也說明了這一點。在垂直分布圖中,混合的貢獻集中在水深200 m以上,并在表層和約50 m處存在兩個峰值。這與海洋上層層結性強是一致的,兩個峰值則可能分別對應于海洋表層由于熱量、淡水輸入形成的強層結和季節性溫躍層形成的強層結。

為進一步分析重力位能改變和層結的關系,本文分析了低緯和高緯兩個樣本浮力頻率和混合致重力位能改變之間的關系(圖3)。低緯樣本取自赤道西太平洋(0.5°S,179.5°E)處,高緯樣本取自南大洋(60.5°S,100.5°W)處。從圖3中可以看 出,無論緯度高低,混合的貢獻和浮力頻率垂直分布形態都有非常好的吻合關系,在高緯二者幾乎重合,低緯則略有出入。原因可能在于垂直混合并不是直接作用于密度,而是分別混合溫鹽進而影響密度,因此海水密度隨溫鹽壓的非線性變化會造成重力位能改變和代表密度層結的浮力頻率不完全一致。

圖2 常混合系數(0.1cm2 s?1)下(a)全球重力位能變化空間分布及(b)混合對重力位能貢獻平均垂直分布

圖3 常混合系數(0.1cm2 s?1)下低緯和高緯混合對重力位能貢獻垂直分布示例

圖4 不同混合系數和浮力頻率下重力位能變化:(a)不同常混合系數下全球大洋重力位能改變;(b)不同混合系數和浮力頻率下的平均重力位能改變

圖4給出了混合系數、浮力頻率、重力位能三者之間的關系。由圖4a可以看出,混合系數和重力位能改變幾乎呈線性關系,在0.1 cms的常混合系數下,全球重力位能改變約為0.08 TW(1TW=10W),而在12 cms的常混合系數下,重力位能改變可達12 TW。假定全球平均混合效率為0.2(Osborn,1980),則等效需要約60 TW的外部能量,這和Wang and Huang(2004a,2004b)的研究認為瞬時風變化能夠輸入海洋約60 TW能量的結論是一致的。密度隨溫鹽壓的非線性效應在這里雖然存在,但并不顯著。這個計算結果和Urakawa and Hasumi(2009)的研究是一致的。由圖4b可知,層結趨于穩定和混合系數的增大都可以增加混合對重力位能的改變,三者之間主要呈正的線性關系,但也存在一定的非線性。結合前面的分析可得出,非線性關系主要存在于層結也即浮力頻率和重力位能改變之間。

4 重力位能改變對風致湍流動能耗散參數化的影響

由上文的分析可知,湍流混合能夠引起大洋重力位能不可忽略的改變,尤其在混合系數為12 cms時,對應所需的外部能量約為60TW,這個量級已經十分接近風(Wang and Huang,2004a)和潮汐(Munk and Wunsch,1998)能夠輸入大洋的能量之和63.5 TW。如此大的能量必然會對大洋中的各種物理過程產生不可忽視的影響,不考慮這種影響的參數化方案可能同實際間存在較大的偏差。

Ledwell et al.(2011)通過DIMS(Diapycnal and Isopycnal Mixing Experiment in the Southern Ocean)試驗得出了在南大洋東南太平洋海盆(62.5°S~55°S,110°W~95°W)處的混合系數和湍流動能耗散率2009年2月至2010年8月的平均垂直分布[見Ledwell et al.(2011)圖1、圖3]。本文用ERA Interim相同區域和時間的風應力和波高數據,結合小節1.2中的方法計算出相應的湍流動能耗散率,并進一步用Osborn(1980)的方案計算出混合系數以對比分析參數化的結果和實際觀測間的差異。計算結果見圖5a。

同Ledwell et al.(2011)的結果對比可以看出,參數化的湍流動能耗散率明顯偏大,在水深100~400 m處偏大2個量級(圖5c),計算出的混合系數也明顯偏大2個量級(圖5b)。本文認為這是沒有考慮到垂直混合會增加重力位能的結果。如果考慮到由風浪輸入表層的湍流動能在向下傳播的過程中會不斷轉化為重力位能,那么參數化的結果應當扣除這部分能量才更合理。根據上文的分析,表層的混合對重力位能的貢獻較小,因此可以假定表層垂直混合不影響混合系數。在此假定下可以計算從表層開始每一層參數化計算出的湍流動能耗散扣除該層以上轉換為重力位能后的結果,從而可得修正后的湍流動能耗散,并進一步計算出該層混合系數。結果見圖5b、c。

由圖5b、c可知,修正后的湍流動能耗散比未經修正直接參數化的結果小2個量級,基本同Ledwell et al.(2011)的觀測相同。但修正后的混合系數仍偏大一個量級,說明在湍流混合層中仍存在未知的能量耗散因素。由圖5c還可看出,垂直混合轉化為重力位能在垂直方向上是一個累積的效應,雖然在風浪攪拌層混合轉化為重力位能的量級遠小于湍流動能耗散,但其累積效應到湍流混合層中已經可以顯著影響到湍流耗散參數化。

圖5 海洋上層修正前后的風浪致湍流動能耗散比較:(a)海洋上層風浪致湍流混合分層示意圖;(b)修正前后的混合系數比較;(c)修正前后的湍流動能耗散率和湍流引起的重力位能改變比較

5 結論

本文利用WOA09資料,統計分析了不同混合系數下的大洋重力位能的變化,并探討了其空間分布及可能原因。在此基礎上,本文進一步討論了垂直混合將湍流能量轉化為重力位能后對湍流參數化的影響。主要的研究結論如下。

(1)垂直混合可以引起顯著的大洋重力位能變化。僅0.1 cms的混合系數就可引起約0.08 TW的重力位能改變,需要約0.4 TW的外部能量,實際大洋中外部能量遠超這個數值,因此在氣候態積分的海洋模式中,混合對重力位能的影響是不可忽視的。重力位能的改變程度同混合系數的大小和層結的強弱有關,一般來說由于海洋中以穩定層結為主,因此混合對重力位能的貢獻以將外部能量轉化為重力位能為主。混合系數越大,層結越強,垂直混合轉化的重力位能越多。

(2)由于密度隨溫鹽壓變化是非線性的,混合增加的重力位能并不和密度層結線性相關,也不嚴格和混合系數大小線性相關。但從全球積分的角度可以認為混合系數和重力位能增加有近似線性關系:。單位為TW,單位為cms。

(3)湍流混合參數化方案應當考慮湍流對外做功的效應。尤其在海洋次表層湍流參數化中,恰當考慮湍流轉換為重力位能效應后(Baumert,2005)的方案可以得到和實際觀測更接近的湍流動能耗散率,進一步計算出的混合系數也比修正前更接近實際觀測。

在大洋中湍流不再僅僅是起耗散作用的能匯,而在一定程度上是能量通道,將外部能量轉換至大洋重力位能。從能量的角度看,大洋中不僅存在傳統意義上的障礙層——溫躍層,也存在能量意義上的障礙層,也即混合轉化重力位能極大區。在能量障礙層,湍流運動因需對外做功而受到抑制,從而進一步影響海洋中物質的湍流混合擴散。提供給混合的外部能量除了因轉化為重力位能消耗外,還可能轉換為其他形式,因此不僅風浪引起的湍流能量耗散參數化應當考慮不同能量間的轉換,混合系數的參數化也應考慮這一點,這些都是值得深入研究的科學問題。此外,考慮能量轉換的參數化方案在數值模式中的具體應用,也是未來值得深入研究的問題。

Baumert H Z. 2005. A novel two-equation turbulence closure for high Reynolds numbers. Part B: Spatially non-uniform conditions [M]// Marine Turbulence: Theories, Observations, and Models. Cambridge: Cambridge Universith Press, 31–43.

Charnock H. 1955. Wind stress on a water surface [J]. Quart. J. Roy. Meteor. Soc., 81 (350): 639–640.

Craig P D, Banner M L. 1994. Modeling wave-enhanced turbulence in the ocean surface layer [J]. J. Phys. Oceanogr., 24 (12): 2546–2559.

Gregg M C. 1989. Scaling turbulent dissipation in the thermocline [J]. J. Geophys. Res., 94 (C7): 9686–9698.

黃瑞新. 1998. 論大洋環流的能量平衡[J]. 大氣科學, 22 (4): 562–574. Huang Ruixin. 1998. On the balance of energy in the ocenic general circulation [J]. Scientia Atmospherica Sinica (in Chinese), 22 (4): 562–574.

Kocsis O, Prandke H, Stips A, Simon A, et al. 1999. Comparison of dissipation of turbulent kinetic energy determined from shear and temperature microstructure [J]. J. Mar. Syst., 21 (1?4): 67–84.

Ledwell J R, Watson A J, Law C S. 1998. Mixing of a tracer in the pycnocline [J]. J. Geophys. Res., 103 (C10): 21499–21529.

Ledwell J R, Montgomery E T, Polzin K L, et al. 2000. Evidence for enhanced mixing over rough topography in the abyssal ocean [J]. Nature, 403: 179–182.

Ledwell J R, St Laurent L C, Girton J B, et al. 2011. Diapycnal mixing in the antarctic circumpolar current [J]. J. Phys. Oceanogr., 41: 241–246.

Mellor G L, Yamada T. 1982. Development of a turbulence closure model for geophysical fluid problems [J]. Rev. Geophys., 20: 851–875.

Munk W, Wunsch C. 1998. Abyssal recipes II: Energetics of tidal and wind mixing [J]. Deep Sea Research Part I: Oceanographic Research Papers, 45: 1977–2010.

Munk W H. 1966. Abyssal recipes [J]. Deep Sea Research and Oceanographic Abstracts, 13: 707–730.

Nagasawa M, Hibiya T, Yokota K, et al. 2007. Microstructure measurements in the mid-depth waters of the North Pacific [J]. Geophys. Res. Lett., 34: L05608.

Osborn T R. 1980. Estimates of the local rate of vertical diffusion from dissipation measurements [J]. J. Phys. Oceanogr., 10: 83–89.

Soloviev A, Lukas R. 2003. Observation of wave-enhanced turbulence in the near-surface layer of the ocean during TOGA COARE [J]. Deep Sea Research Part I: Oceanographic Research Papers, 50: 371–395.

St Laurent L, Naveira Garabato A C, Ledwell J R, et al. 2012. Turbulence and diapycnal mixing in Drake Passage [J]. J. Phys. Oceanogr., 42: 2143–2152.

Terray E A, Donelan M A, Agrawal Y C, et al. 1996. Estimates of kinetic energy dissipation under breaking waves [J]. Journal of Physical Oceanography, 26: 792–807.

Thompson A F, Gille S T, MacKinnon J A, et al. 2007. Spatial and temporal patterns of small-scale mixing in Drake Passage [J]. J. Phys. Oceanogr., 37: 572–592.

Urakawa L S, Hasumi H. 2009. The energetics of global thermohaline circulation and its wind enhancement [J]. J. Phys. Oceanogr., 39: 1715–1728.

Wang W, Huang R X. 2004a. Wind energy input to the surface waves [J]. J. Phys. Oceanogr., 34: 1276–1280.

Wang W, Huang R X. 2004b. Wind energy input to the Ekman layer [J]. J. Phys. Oceanogr., 34: 1267–1275.

張宇, 林一驊, 王輝. 2014. 垂直湍流輸送對大洋的重力位能和混合過程的影響 [J]. 大氣科學, 38 (5): 838?844, doi:10.3878/j.issn.1006-9895. 1401.13105. Zhang Yu, Lin Yihua, Wang Hui, et al. 2014. Impact of vertical turbulence on ocean gravitational potential energy and the tracer mixing process [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 38 (5): 838?848.

Impact of Vertical Turbulence on Ocean Gravitational Potential Energy and the Tracer Mixing Process

ZHANG Yu, LIN Yihua, and WANG Hui

1,100081;2,,,100029

Using a WOA09 (World Ocean Atlas 2009) data set of objectively analyzed in situ temperature and salinity, we calculate ocean gravitational potential energy (GPE) and investigate the relationships between the mixing coefficient, buoyancy frequency, and GPE. On that basis, we further explore the impact of energy conversion, caused by turbulent mixing, on turbulent parameterization. The research shows that ocean vertical turbulence is not only a kinetic energy sink but also an important way of external energy transformation to GPE. When the mixing coefficient is 0.1 cms, GPE will increase 0.08 TW. Based on the results of other authors, we conclude that external energy can induce a global average mixing coefficient of up to 12 cms. In general, the more stable the stratification of the ocean and the larger the mixing coefficient, the higher the GPE increase. Parameterization can obtain more realistic results with careful treatment of turbulent kinetic energy conversion to GPE.

Mixing coefficient, Gravitational potential energy, Turbulent kinetic energy dissipation rate, Turbulent mixing parameterization

1006?9895(2014)05?0838?07

P731

A

10.3878/j.issn.1006-9895.1401.13105

2013?01?07,2014?01?23收修定稿

國家自然科學基金資助項目41175058

張宇,男,1984年出生,博士研究生,主要從事大氣邊界層和大洋環流數值模擬研究。E-mail: zhangyu7305@gmail.com

林一驊,E-mail: linyh@lasg.iap.ac.cn

主站蜘蛛池模板: 亚洲色无码专线精品观看| 免费无码在线观看| 久久亚洲美女精品国产精品| 福利视频久久| 午夜福利在线观看入口| 国产情侣一区| 在线看免费无码av天堂的| 黄色网页在线播放| 国产第一页屁屁影院| 国产精品久久久久久久久久98| 91成人在线观看| 在线免费不卡视频| 91成人在线观看| 久久国产精品夜色| 亚洲AV无码乱码在线观看裸奔 | 久久精品亚洲热综合一区二区| 国产黄色免费看| 亚洲精品中文字幕午夜| 影音先锋亚洲无码| 日韩精品免费一线在线观看| 91小视频在线播放| 亚洲视频四区| 伊人成人在线| 精品無碼一區在線觀看 | 91丝袜在线观看| 啪啪免费视频一区二区| 国产在线自在拍91精品黑人| 精品无码日韩国产不卡av| 久久亚洲日本不卡一区二区| 巨熟乳波霸若妻中文观看免费| 九色视频线上播放| 免费无码一区二区| 久久精品亚洲中文字幕乱码| 综合色婷婷| 欧美日韩在线亚洲国产人| 97综合久久| 国产一在线| 亚洲美女操| 亚洲无卡视频| 人妻丰满熟妇αv无码| 亚洲精品免费网站| a亚洲天堂| 免费精品一区二区h| 国产粉嫩粉嫩的18在线播放91| 精品丝袜美腿国产一区| 色偷偷男人的天堂亚洲av| 久久精品午夜视频| 免费jjzz在在线播放国产| 久久久久无码国产精品不卡| 青青网在线国产| 四虎国产在线观看| 欧美特黄一级大黄录像| 亚洲区欧美区| 五月激情婷婷综合| 色视频国产| 无码电影在线观看| 91福利一区二区三区| 国模视频一区二区| 动漫精品啪啪一区二区三区| 久久99精品久久久久纯品| 亚洲自拍另类| 色AV色 综合网站| 亚洲第一区在线| 久久精品丝袜高跟鞋| 国产剧情无码视频在线观看| 国产成人高清在线精品| 五月天综合婷婷| 精品人妻一区无码视频| 99ri国产在线| 国产精品片在线观看手机版| 国产精品不卡片视频免费观看| 在线亚洲精品自拍| 国产精品亚欧美一区二区| 国产午夜无码专区喷水| 国产欧美精品一区二区| www.亚洲国产| 91色国产在线| 在线视频亚洲色图| 青青青伊人色综合久久| 免费在线视频a| 青草视频网站在线观看| 亚洲欧美不卡中文字幕|