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

中緯度地磁暴期間熱層垂直風響應機制的模擬*

2022-04-13 03:23:00蘇燁李婧媛呂建永王明魏官純孫夢熊世平李正
空間科學學報 2022年2期
關鍵詞:水平

蘇燁 李婧媛 呂建永 王明 魏官純 孫夢 熊世平 李正

(南京信息工程大學大氣物理學院空間天氣研究所 南京 210044)

0 引言

熱層大氣位于中間層與逸散層之間,是中高層大氣的重要組成部分,熱層大氣風場受地磁活動影響很大。相對于水平風而言,熱層垂直風很小,測量的難度也導致相關觀測結果非常缺乏[1]。但是對地磁暴期間垂直風場的研究從未停止,關于水平風場與垂直風場之間聯系的理論研究可以追溯到20 世紀70年代,Dickinso等[2]認為連續性方程中的垂直速度可分為由壓力面總質量散度決定的垂直速度和等壓面抬升引起的垂直速度。Rishbeth等[3]將總的垂直風速也分為兩部分:表示在高度上平衡水平輻合或輻散的散度速度;表示熱層收縮或膨脹時等壓面上升或下降的氣壓速度。Burnside等[4]建立了關于垂直風場與水平風場散度相關方程。由于高層大氣探測手段缺乏,對于垂直風場與水平風場之間聯系的研究僅停留在理論階段。

隨著天基和地基探測的發展,關于地磁暴期間熱層垂直風場變化的分析也隨之增多,其中地基探測數據多使用法布里–珀羅探測儀(FPIs)研究熱層的垂直風場。Hernandez[5]指出在準穩態或隨時間變化的情況下,強烈的局部加熱和地磁擾動造成的熱層風場加速可產生接近100 m·s–1以上的垂直風速。Spencer等[6]使用WATS(Wind and Temperature Spectrometer)觀測數據得到,在有地磁擾動的情況下,極光卵附近300 km 以上的熱層區域能夠觀測到的最大垂直風速約100~200 m·s–1。Harding等[7]分析了北美熱層–電離層觀測網(NATION)5年(2011-2016年)的數據,表明在地磁暴期間中緯度地區能觀測到風速較大的熱層垂直風變化,并且這種垂直風的變化會持續存在。Hu等[8]分析了從北極黃河站全天FPIs 觀測中獲取的熱層垂直風數據,得到高熱層與低熱層的垂直風平均速度分別為40 m·s–1和15 m·s–1,并且存在日變化和與地磁活動水平相關的變化。以上研究均證明在地磁暴期間,由地磁擾動引起的熱層垂直風場變化非常劇烈,但是造成磁暴期間垂直風場變化的物理機制尚不明了。

FPIs 數據和非相干散射雷達(ISR)數據的分析研究印證了垂直風場與水平風場的聯系[9,10]。Biondi[11]研究表明,觀測到的向上和向下垂直運動是對中緯度熱層水平風場輻散和輻合的響應。Rees等[12]發現高度積分(height-integrated)的垂直風速可能超過100~150 m·s–1,垂直風的變化與水平風場的高度相關,這種變化是對地磁能量輸入的響應,平均垂直風的強度能夠達到約30 m·s–1(或100 km·h–1)且持續3~6 h,認為這種持續且強烈的垂直風場變化是熱層受地磁強烈擾動的表現。Peteherych等[13]的研究結果表明,在110~125 km 高度,焦耳加熱直接導致上行風發生,而水平風的散度導致了下行風。Smith等[14]發現在地磁擾動強且伴隨強烈西風時,觀測到高達150 m·s–1的向下氣流,而夜間這種垂直氣流則以相同速度上行。他們認為造成這種較高熱層垂直風場變化的一個主要機制是水平風的散度。由此可以得出,在地磁暴期間,熱層中緯度垂直風的變化與水平風密切相關。

近年也有研究利用模型研究磁暴期間垂直風場變化與溫度變化的聯系,Li等[15]使用TIMEGCM 模型研究了中緯度中間層-低熱層(MLT)中性溫度對地磁暴的響應。結果表明地磁暴發生后,在中緯度MLT區域的垂直風先發生變化,這種垂直風的變化通過絕熱加熱/冷卻與垂直熱輸送過程導致局地暴時溫度變化。Li等[16]進一步研究發現,中緯度地區來自熱層的下行垂直風會造成MLT 區域的溫度增加,這種溫度增加在磁暴早期會驅動局地氣流向赤道運動,從而造成向南的擾動風場。隨著磁暴的發展,MLT 溫度升高會在中緯度產生分流,使得擾動風的方向變為向北向上。一直以來,由于探測手段和模型精度的制約,對于磁暴期間垂直風場變化的研究十分缺乏。雖然垂直風在量級上比水平風小得多,但對于熱層大氣環境(溫度及環流等)的影響卻是顯著的,在垂直方向梯度較小的波動也能引起熱層巨大的變化,因此研究磁暴期間垂直風場變化對于理解太陽活動對中高層大 氣影響具有重要意義。

1 數據和方法

1.1 TIMEGCM 模型

研究主要基于TIMEGCM 模型。TIMEGCM模型的發展可追溯至20 世紀80年代初,美國國家大氣研究中心(NCAR)開始建立TGCM 模型(Thermospheric General Circulation Models),以研究熱層的全球溫度、環流和化學成分及其對太陽和極光活動的響應,隨后擴展為三維TIGCM 模型,包含電離層與熱層之間電動力學相互作用自洽的TIEGCM 模型。TIEGCM 模型向下延伸至30 km 高度,包括中間層和平流層上層,被稱為TIMEGCM,其需要的輸入參數包括57 個主要波長的太陽紫外通量(由F10.7通量參數化得到)、極光粒子沉降、高緯度電場以及從30 km 下邊界向上傳播的潮汐,TIMEGCM 在了解高層大氣的特征方面發揮了重要作用[17-22]。在本研究中,使用的TIMEGCM 模型水平分辨率為2.5°×2.5°(緯度×經度),垂直分辨率為1/4 個標高。高緯對流電場使用由3 h的Kp指數驅動的高緯對流電場 Heelis 模型[23]。

1.2 垂直風診斷

中性大氣連續性方程為

在Boussinesq 假設條件下可得

化簡得到

其中,un和vn分別為緯向和經向風速,R為地球半徑,λ為緯度,φ為經度,z為高度,W為垂直風速。等價于,表示緯向風場沿著模型水平網格的變化等價于,表示經向風場沿模型水平網格的變化;表示垂直風沿模型高度網格的變化。因此可以得到

通過式(4)可以看出,垂直風速由高層的垂直風速變化及同層水平風速變化求得。由于TIMEGCM模式計算的是壓力面間的垂直風速,因此地面加熱使得壓力面起伏產生的垂直風速不包括在式(4)的計算內。

使用TIMEGCM 模型進行敏感性實驗:一個由地磁平靜條件的2.3 常數Kp值(圖1 中黑色橫虛線)驅動,另一個由真實的Kp值(圖1 中黑色折線)驅動。3 h 指數Kp值在2005年9月10日為2+、3?、4、3+、4、6?、5+、5。可以看到,在24 h的Kp值變化中,15~18 h的Kp值最大。

圖1 2005年9月10日全天的Kp 值(黑色折線為Kp 值,黑色橫虛線表示Kp=2.3,黑色豎虛線表示磁暴擾動開始的時間)Fig.1 Kp value of the whole day on 10 September 2005 (The black polyline is the Kp value,the black horizontal dashed line represents Kp=2.3,and the black vertical dashed line is the time of the start of the magnetic storm)

除了Kp值條件不同外,這兩次實驗在00:00 UT使用完全相同的初始條件和相同的其他模型驅動參數。因此,兩次實驗結果之間的差值即為地磁暴對風場的影響。其中,06:00 UT(黑色豎虛線)為2005年9月10日地磁暴擾動的開始時間[24]。

2 結果與分析

以下涉及的風場及溫度數據均采用磁暴值減去平靜值的方法處理,確保是地磁暴期間由于地磁擾動影 響所引起的變化。

2.1 水平風場沿水平網格的變化與垂直風場在高度上的對比

圖2 顯示TIMEGCM 模擬了中緯度(30°N-50°N 緯度平均)地磁暴期間,4 個時間點(06:00 UT、10:00 UT、14:00 UT、18:00 UT)的熱層(100~650 km高度)水平風場沿經度和壓力面網格的變化(圖2 a~d)和垂直風場(圖2 e~h)。可以看出,06:00 UT的風場整體速度較小,即由磁場擾動引起的風速變化較小,且水平風場與垂直風場的變化基本均為負值,此時主要是較高高度(約250 km 以上)的垂直風拉動較低高度(約250 km 以下)的垂直風且保持向下傳播;而到了10:00 UT,磁場擾動增大,風場的整體擾動速度隨之增大,出現了向上的垂直風,但是擾動的垂直風場始終保持其連續性;到了14:00 UT 和18:00 UT,磁場擾動劇烈,風場速度變化激增,但依然是較高高度拉動較低高度。并且在這4 個圖中,250 km 高度以上的水平風場散度與垂直風場始終保持基本一致。

圖2 中緯度地磁暴期間熱層水平風場沿緯度壓力面網格的變化(F)與垂直風場(向上為正,向下為負;黑色虛線表示250 km 高度,黃色和綠色虛線分別代表06:00 UT 和18:00 UT 時刻)Fig.2 Variation of horizontal wind along the grid of longitude (F) and vertical wind in the thermosphere during mid-latitude geomagnetic storms (Up is positive,down is negative;the black dotted line indicates 250 km height,and the yellow and green dashed lines indicate 06:00 UT and 18:00 UT respectively)

圖3 顯示TIMEGCM 模擬了中緯度(30°N-50°N)地磁暴期間,4 個時間點(06:00 UT、10:00 UT、14:00 UT、18:00 UT)的熱層(100~650 km)水平風場沿緯度和壓力面網格的變化(圖3 a~d)和垂直風場(圖3 e~h)在107.5°W。可以看出,在4 個時間段內,250 km 高度以上水平風場的散度與垂直風場始終保持基本一致。

圖4 顯示TIMEGCM 模擬了中緯度(30°N-50°N)地磁暴期間,4 個時間點(06:00 UT、10:00 UT、14:00 UT、18:00 UT)的熱層(100~650 km 高度)水平風場在92.5°E 沿緯度和壓力面網格的變化(圖4 a~d)和垂直風場(圖4 e~h)。可以看出,在4 個時間段內,250 km 高度以上水平風場的散度與垂直風場始終保持基本一致。由圖3 和圖4 可得,固定經度下的緯度變化并不影響本節的結論。

圖3 中緯度地磁暴期間熱層水平風場沿緯度壓力面網格的變化(F)與垂直風場 (向上為正,向下為負;黑色虛線為250 km)Fig.3 Variation of horizontal wind along the grid of longitude (F) and vertical wind in the thermosphere during mid-latitude geomagnetic storms (up is positive,down is negative,and the black dotted line is marked as 250 km)

圖4 中緯度地磁暴期間熱層水平風場沿緯度壓力面網格的變化(F)與垂直風場 (向上為正,向下為負;黑色虛線為250 km 高度)Fig.4 Variation of horizontal wind and vertical wind along the grid of longitude (F) in the thermosphere during mid-latitude geomagnetic storms(Up is positive,down is negative,and the black dotted line is marked as 250 km height)

結果表明,在250 km 高度以上,水平風場的散度與垂直風場基本相同,即可以使用水平風場散度在地磁暴期間的變化解釋垂直風場變化。在250 km 高度以下,可以看到垂直風場的擾動變化都是由于較高高度擾動垂直風場拉動的,低高度垂直風場的擾動變化量級遠小于高高度,因此100~650 km 高度的擾動垂直風場具有很好的連續性。由此可以進一步研究是經向還是緯向風場在地磁暴期間的變化導致了垂直風 場的產生。

2.2 經向與緯向風場對水平風場變化的影響

將水平風場拆分成緯向風場與經向風場,將三個風場在固定時間點及統一坐標軸的情況下進行比較,通過數據對比可以得到兩個風場對水平風場變化的貢獻,進一步理解影響垂直風場變化的物理機制。

圖5 顯示TIMEGCM 模擬了中緯度(30°N-50°N)地磁暴期間(高度平均約250 km 以上)熱層緯向風場散度(圖5 a~d)、經向風場散度(圖5 e~h)及水平風場的散度(圖5 i~l)沿經緯度網格在4 個時間點(06:00 UT、10:00 UT、14:00 UT、18:00 UT)上的變化。地磁暴擾動開始時(06:00 UT),向東變化的緯向風與向南變化的經向風疊加合成后產生了向下的垂直風,這時明顯可以看出經向風的變化起到主導作用;10:00 UT 時,垂直風主要集中在夜側(06:00 UT 以西18:00 UT 以東的區域),變化更大的緯向風與經向風疊加合成后,垂直風的變化主要由緯向風的散度造成;到了14:00 UT,其與10:00 UT 的疊加合成效應相同,可以看出緯向風的變化對垂直風產生了更大貢獻;在18:00 UT,通過分析夜側風場,可以看出緯向風變化對垂直風場變化的影響與經向風相差不大,但在日側(06:00 UT 以東18:00 UT 以西的區域),主要是緯向風變化影響了垂直風。

圖5 中緯度地磁暴期間高熱層緯向風場、經向風場與水平風場的散度暴時變化(向東和向北為正,向西和向南為負;黃色虛線和綠色虛線分別注明了06:00 UT 和18:00 UT)Fig.5 Divergence velocity changes of the zonal wind,meridional wind and horizontal wind during geomagnetic storm in the middle latitude (the east is positive and the west is negative;the yellow dotted line and the green dotted lines indicate 06:00 UT and 18:00 UT respectively)

因此,擾動開始時,經向風變化比緯向風場更大一些,即對于水平風場散度的影響也更大。隨著地磁暴的發展,在10:00 UT 和14:00 UT,緯向風變化相對于經向風場更大一些,同理可得,在這個時間對于水平風場散度影響更大的是緯向風場。到了18:00 UT,地磁擾動更為劇烈,緯向風場與經向風場的變化也隨之加劇。隨著兩個風場進一步的疊加和合成作用,水平風場的變化也同樣變得更為劇烈,但依舊是緯向風變化相對較大,此時其對水平風場散度變化的作用也更大。因此可以得到如下結論:在地磁暴擾動開始階段(06:00-08:00 UT),在250 km 高度以上,相比緯向風場,熱層經向風場的變化是產生垂直風的更大動力源,而當地磁暴擾動發展起來后,緯向風場的變化相 對更大,在垂直風的變化趨勢上起到更大推動作用。

2.3 水平風場與溫度傳播的聯系

為了將經向風和緯向風沿水平方向的變化與其風場本身及溫度的變化聯系起來,將經向風、緯向風及溫度沿經緯線的變化與風場自身數據在固定時間點及統一坐標軸的情況下進行比較,由此得到地磁暴期間溫度傳播與水平風場變化的內在聯系。

圖6 顯示TIMEGCM 模擬了中緯度(30°N-50°N)地磁暴期間(高度平均250 km 以上)熱層溫度場(圖6 a~d)和緯向風場(圖6 e~h)沿緯線的差值以及緯向風場(圖6 i~l)在4 個時間點(06:00 UT、10:00 UT、14:00 UT、18:00 UT)上的變化。在06:00 UT,溫度梯度(東邊格點溫度減去西邊相鄰格點溫度)值極小,且主要出現在夜側;10:00 UT 夜側的溫度梯度正值變化區域比負值的大一些,日側也出現了溫度梯度負值區域;在14:00 UT 溫度梯度值比之前增大一些,并且在日側主要為正值區域,在夜側主要為負值區域;在18:00 UT,相比之前的3 個時間點,溫度梯度值最大,延續了14:00 UT 的正負值區域的分布特點。圖6 中?u的變化與圖5中的變化基本一致。緯向風場在06:00 UT 主要為夜側的西向風,但在其東側也有一個較小的東向風;10:00 UT 夜側西向風增大,東側的東向風也隨之增大,并且在日側也出現了明顯的西向風,但量級比夜側的風場小得多;到了14:00 UT,100°W 東側(日側)和西側(夜側)均存在一個較小的東向風,其余區域均為風速較大的西向風,且西向風的風速大值區域主要存在于夜側;18:00 UT 整個中緯度全是風速極大的西向風,但夜側的西向風風速大于日側。

圖6 中緯度地磁暴期間熱層緯向風場和溫度場沿緯線的差值與緯向風場(向東為正,向西為負;黃色和綠色虛線分別表示06:00 UT 和18:00 UT 時刻)Fig.6 Difference between the zonal wind and temperature along the latitude line and the zonal wind during geomagnetic storm in the middle latitude (The east is positive,the west is negative,and the yellow and the green dotted lines indicate 06:00 UT and 18:00 UT respectively)

地磁擾動開始(06:00 UT)后,溫度梯度值在中緯度夜側比日側大得多,沿東向緯線升高的溫度(圖6 a~d)產生了西向的緯向風(圖6 i~l),反之沿東向緯線降低的溫度產生東向的緯向風。到了10:00 UT,溫度沿東向緯線在夜側的變化比日側大得多,仍然是向東增加的溫度梯度形成西向風,向西增加的溫度梯度形成東向風。隨著擾動增強,到了14:00 UT,日側溫度梯度也隨之增強,但依然是夜側溫度梯度更大,緯向風場仍遵循與溫度梯度的關系。進入18:00 UT,日側和夜側的緯向風場整體速度都很大,但隨著擾動進一步增強,溫度梯度與緯向風場依舊呈現出很好的對應性,溫度梯度在日側呈現出沿東向緯線升高的趨勢,在夜側呈現出沿東向緯線降低的趨勢,而緯向風場無論在日側還是夜側都是西向風。可以看出,溫度梯度和緯向風場均隨地磁暴擾動的增強而增大。并且在地磁暴擾動期間,溫度梯度與緯向風的變化趨勢基本一致,可以看出沿緯線分布的緯向風基本由沿緯線分布的溫度梯度決定。

圖7 顯示TIMEGCM 模擬了中緯度(30°N-50°N)地磁暴期間(高度平均250 km 以上)熱層溫度場(圖7 a~d)和經向風場(圖7 e~h)沿經線的差值以及經向風場(圖7 i~l)在4 個時間點(06:00 UT、10:00 UT、14:00 UT、18:00 UT)上的變化。溫度梯度(北南)在06:00 UT 的值極小,且主要是沿北向經線增大的溫度梯度(正)出現在夜側;10:00 UT 夜側的溫度梯度正值區域擴大,在日側出現一個較小的沿北向經線減小的溫度梯度(負)區域;14:00 UT 的溫度梯度正值區域急劇增大并擴展到整個北半球,但夜側的溫度梯度量級更大;18:00 UT 的溫度梯度正值在整個中緯度的量級急劇增大,此時日側和夜側的量級相差不大。圖7 中?v的變化與圖5中的變化基本一致。經向風場在06:00 UT 僅在夜側產生一個較小的南向風;10:00 UT 的夜側南向風風速有所增大;14:00 UT 的經向風場延續了上一個時間點的夜側南向風,同時在日側也出現一個風速較小的南向風;南向風在18:00 UT 擴展至整個北半球,并且風速急劇增大,形成一個超大型的南向風。

圖7 中緯度地磁暴期間熱層溫度場和經向風場沿經線的差值與經向風場(向東為正,向西為負,黑色虛線標注的是250 km,黃色和綠色虛線分別表示06:00 UT 和18:00 UT 時刻)Fig.7 Difference between the temperature of the hyperthermia and the meridional wind along the meridian during the mid-latitude geomagnetic storm and the meridional wind (North is positive,south is negative,and the yellow and the green dotted lines indicate 06:00 UT and 18:00 UT respectively)

地磁擾動開始(06:00 UT)后,溫度梯度在中緯度夜側的量級大于日側,沿北向經線升高的溫度(圖7a~d)產生了南向的經向風(圖7 i~l),即高緯度的溫度沿南向的經向風向低緯度擴散,且僅有南向的經向風。在10:00 UT,溫度梯度沿北向經線在夜側的變化占溫度傳播的主導地位,依然遵循沿北向經線升高的溫度產生南向經向風的規律。隨著擾動增強,到了14:00 UT,日側溫度梯度沿北向經線的變化也隨之增強,但依然是夜側的溫度梯度更大,經向風場仍遵循與溫度梯度的關系,整個中緯度北半球都是南向經向風。而進入18:00 UT,日側和夜側的經向風場整體速度都很大,但隨著擾動進一步增強,溫度梯度的變化與經向風場對應得更好,溫度梯度在整個中緯度北半球呈現出沿北向經線急劇升高的趨勢,即越來越多的高緯度溫度沿經線向低緯度傳播,而經向風場都是南向的經向風。因此可以得出,地磁擾動開始后,溫度開始沿著經線擴散,溫度與經向風的變化也隨擾動增大而增強,沿北向經線升高的溫度梯度產生了向南的經向風。因此,溫度梯度對于經向風場的影響在擾動增強時更明顯。

由圖6 和圖7 可知,前兩個小時經向風的變化更明顯地導致了垂直風的變化,由磁暴導致的增溫從高緯度向低緯度擴散,沿經線的風速變化較大;08:00 UT之后,溫度沿緯線擴散得更快,風速沿緯線變化更大,因此緯向風的變化更明顯導致了垂直風的變化。而隨著擾動增強,溫度不再是影響水平風場的唯一因素,緯向風沿緯線的變化比經向風沿經線的變化更大,緯向風場的變化更明顯地導致了垂直風場的變 化。

3 結論

利用TIMEGCM 模型,研究了2005年9月10日中緯度地磁暴期間熱層(100~650 km 高度)水平風場變化對垂直風的影響,依據模型中垂直風求解公式,對數據進行了處理與分析,主要得到結果和結論如下。

(1)在地磁暴期間,250 km 以上的垂直風場變化與水平風場的散度相同,因此可以得到水平風場變化引起垂直風場變化的結論,由于100~650 km 高度的垂直風是連續的,250 km 以下的垂直風場主要由較高高度垂直風拉動。

(2)在地磁暴初期(06:00-08:00 UT),經向風場的變化比緯向風場更大,因此對垂直風場變化做出了更大貢獻,而隨著擾動的增強,緯向風場變化更大,在08:00 UT 之后,相比經向風場,緯向風場對垂直風場的變化影響更大。

(3)結合溫度沿經緯線的變化分析,在地磁暴初期(06:00-08:00 UT),溫度沿經線傳播得更快,經向風沿經線的方向變化更大,而在08:00 UT 之后,溫度沿緯線傳播得更快,緯向風沿緯線的方向變化更大,與結論2 有很好的對應。

雖然已通過TIMEGCM 模型模擬研究了中緯度地磁暴期間熱層垂直風場的變化,以及其與水平風場之間的內在聯系,并得出了初步結論,但是對于地磁暴期間影響熱層水平風場變化的物理過程尚不了解,今后將對其進行更深入研究。

致謝地磁活動指數Kp數據來源于ftp://ftp.ngdc.noaa.gov/STP/GEOMAGNETIC_DATA/INDICES/KP_AP/網站。

猜你喜歡
水平
張水平作品
作家葛水平
火花(2019年12期)2019-12-26 01:00:28
深化精神文明創建 提升人大工作水平
人大建設(2019年6期)2019-10-08 08:55:48
加強上下聯動 提升人大履職水平
人大建設(2019年12期)2019-05-21 02:55:32
水平有限
雜文月刊(2018年21期)2019-01-05 05:55:28
加強自身建設 提升人大履職水平
人大建設(2017年6期)2017-09-26 11:50:44
老虎獻臀
中俄經貿合作再上新水平的戰略思考
建機制 抓落實 上水平
中國火炬(2010年12期)2010-07-25 13:26:22
做到三到位 提升新水平
中國火炬(2010年8期)2010-07-25 11:34:30
主站蜘蛛池模板: 亚洲中文字幕97久久精品少妇| 日本黄网在线观看| 亚洲精品天堂自在久久77| 国产视频a| 色欲色欲久久综合网| 国产拍揄自揄精品视频网站| 亚洲天堂视频在线播放| 亚洲精品国产成人7777| 久草国产在线观看| 一级做a爰片久久毛片毛片| 一级毛片在线直接观看| 青青青国产在线播放| 国产成人夜色91| 一边摸一边做爽的视频17国产| 熟女日韩精品2区| 激情無極限的亚洲一区免费| 色哟哟色院91精品网站 | 久久综合AV免费观看| 国产人成在线视频| 国产在线视频欧美亚综合| 精品欧美一区二区三区久久久| 日韩国产一区二区三区无码| 99人体免费视频| 亚洲成人在线网| 国产另类乱子伦精品免费女| 在线观看热码亚洲av每日更新| 亚洲日韩第九十九页| 国产亚洲欧美日韩在线一区二区三区 | 中文字幕乱码二三区免费| 亚洲欧美一区二区三区图片| 免费啪啪网址| 国产精品无码翘臀在线看纯欲| 国产精品真实对白精彩久久 | 欧美一级黄色影院| 国产成人91精品免费网址在线| 国产欧美在线观看精品一区污| 国产在线观看一区精品| 国产一级在线观看www色| 青草视频网站在线观看| 国产白丝av| 香蕉在线视频网站| 精品国产自在现线看久久| 蜜桃臀无码内射一区二区三区| 欧美一区二区人人喊爽| 国产香蕉国产精品偷在线观看| 国产精品浪潮Av| 五月婷婷综合色| 人妻少妇乱子伦精品无码专区毛片| 激情综合激情| 丁香婷婷久久| 国产无码精品在线播放| 国产综合精品日本亚洲777| 国产成人调教在线视频| 国产一二视频| 国产成人免费| 亚洲人成高清| 久久久久国产精品免费免费不卡| 91口爆吞精国产对白第三集| 国产99在线| 少妇人妻无码首页| 亚洲综合日韩精品| 亚洲最大福利视频网| 国产精品永久免费嫩草研究院| 国产精品一区在线麻豆| 伊人色在线视频| 极品国产一区二区三区| 一本大道无码日韩精品影视| 免费jizz在线播放| 国产日产欧美精品| 国产精品视频猛进猛出| 国产成人精品免费视频大全五级| 久久9966精品国产免费| 国产超碰在线观看| 99久久国产综合精品2020| 亚洲人成影视在线观看| 亚瑟天堂久久一区二区影院| 欧美午夜网站| 免费无码网站| 激情六月丁香婷婷四房播| 国产99久久亚洲综合精品西瓜tv| 又爽又黄又无遮挡网站| 国产精品视频久|