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

中小尺度對流系統的高分辨率 數值模擬近況和未來挑戰

2021-11-06 10:10:34黎慧琦張大林

黎慧琦 張大林,

(1 中國氣象局廣州熱帶海洋氣象研究所,廣州510640;2 中國氣象科學研究院災害天氣國家重點實驗室,北京100081;3 美國馬里蘭大學大氣海洋系,馬里蘭 20742)

0 引言

自從1922年Richardson提出以流體力學為基礎的數值天氣預報概念,數值模式逐漸成為大氣科學中一個重要的領域。尤其20世紀80年代以來,隨著高性能計算機的迅速發展,模式分辨率不斷提高,模式次網格尺度物理過程的逐步改進,加上越來越多觀測資料的同化,大大提高了全球和區域數值天氣預報精度,使得數值天氣預報已成為現代天氣預報業務的基礎和天氣預報業務發展的主流方向。可以說,數值天氣預報所取得的成就是近百年來物理學科各領域中最具革命性的。此外,這些成就還帶動了海洋、水文、環境、生態預測或預報等其他領域的發展。同時,數值模式作為一項重要的工具,能為天氣現象機理研究提供動力學、熱力學一致的四維數據集,以補充觀測網未能捕捉的信息。借助數值模式還能開展敏感性試驗和動力、熱力診斷分析,區分不同的物理過程或影響因子的作用。

業務模式的數值天氣預報和物理機制研究為主所進行的數值模擬是有差別的。對于后者又分為實際個例模擬和理想數值試驗。實際個例模擬采用實際的初始和側邊界條件,通常是從全球再分析數據獲得,模擬的結果需要先和觀測進行對比驗證,然后用于特定問題的研究。理想數值試驗使用隨高度變化的標準探空啟動試驗。而業務模式和數值模擬研究的主要目的也不同。前者主要用于預報氣溫、濕度、風和降水等氣象要素,后者主要研究天氣系統的三維結構、演變及其物理機制。總之,深入理解天氣演變機理和提高氣象預報(模擬)能力是發展數值天氣預報模式的兩個主要目標。

到20世紀末,對積云對流、云微物理、邊界層過程等相關物理過程認識有所加深,相應的參數化方案得以不斷完善,同時得益于計算機的高速發展,數值模式的分辨率已得到顯著提高。早期的業務模式分辨率在380 km左右,而到90年代,區域模式已達5~10 km。眾多業務中心逐漸從采用原始方程模式發展到采用非靜力模式。21世紀以來,不少國家地區的非靜力業務模式的水平分辨率已逐步達到1~3 km,例如美國使用3 km的高分辨率快速更新系統(HRRR),英國使用1.5 km版本的MetUM,瑞士使用1.1 km的COSMO模式,我國GRAPES-Meso模式分辨率也發展至3 km。

近十年來,越來越多的研究嘗試使用分辨率小于1 km,甚至幾十米的大渦數值模擬來探討更精細的物理過程,這些研究所用的數值模式包含的可供選擇的物理方案比業務模式豐富。本文將對中、小尺度對流系統(如颮線、中尺度對流復合體、熱帶氣旋和鋒面系統中的對流、龍卷風暴、超級單體)的高分辨率(主要關注1 km左右以及次千米級)模擬涉及的重要方面進行梳理。第1~3節分別討論在對某些中、小尺度對流系統進行數值模擬時,模式分辨率從12 km提高到1 km以下的作用,改善初始條件的方法及重要性,以及不同分辨率下使用不同物理方案的影響。第4節給出幾個1 km和次千米級的高分辨率模擬例子。最后一節總結中、小尺度對流系統高分辨率數值模擬面臨的挑戰和未來展望。

1 模式分辨率的影響

大氣運動從上萬千米行星環流到毫米量級湍流輸送,橫跨12個量級尺度。隨著計算機能力的不斷提升,人們總是設法模擬更小尺度的大氣現象,這就使得數值模擬的分辨率不斷提高。一般來說,要模擬更小尺度的天氣現象或梯度大的氣象變量,尤其由復雜地形、非均勻地表所產生的大氣運動,需要使用更高的分辨率。我們可以認為,這主要因為使用高分辨率網格距雖然仍存在次尺度物理參數化的問題,但可減小由重要物理過程(如對流、邊界層)參數化產生的嚴重不確定性,從而更好地顯示可分辨尺度各物理定律對大氣運動的控制度。在當前的計算機條件下,大多高分辨率數值模擬仍需要采用多重嵌套網格。在使用同樣的初始條件和物理方案的前提下,主要有以下兩個途徑改善業務預報:一個是通過提高分辨率改善中、小尺度對流系統的確定性預報結果,另一個是采用更多成員的集合預報以避免單個模式預報的不確定性。業務預報單位往往需要考慮重點將計算資源分配在哪個途徑。而對數值模擬研究來說,大量研究顯示提高水平分辨率有利于改善地形環流,城、郊以及水、陸差異,中尺度對流系統的結構,臺風和降水的強度等的模擬。可以認為,提高分辨率的正效果也可能與復雜微物理方案模擬得到網格尺度的潛熱釋放有關。

Weisman等采用12~1 km的分辨率模擬颮線,發現分辨率越高,模擬的颮線結構和演變更接近真實。Bryan等分析使用更高分辨率(1~0.125 km)對颮線模擬的影響,他們的研究發現1 km的模擬不能合理重現湍流動能特征,而更高分辨率的模擬(如0.125 km)能更好地描述湍流特征(圖1),這使得使用兩種分辨率所模擬的颮線在系統相速度、 云層厚度、對流組織化、降水量等方面有明顯差異。Bryan等利用4、1、0.25 km三種分辨率對颮線的模擬研究指出,較低分辨率的模擬由于模擬的對流單體偏大,缺乏夾卷中層空氣,使得云頂發展得偏高,產生偏多降水,而較高分辨率的模擬由于可分辨湍流的發展,有更多的云水蒸發。在很多情況下,分辨率主要是影響上升氣流的模擬,從而影響水汽凝結或云水蒸發過程。Verrelle等使用中層較濕的初始探空啟動模擬,則發現總降水量隨著分辨率的降低而減少。Potvin等對超級單體的模擬研究顯示4 km分辨率不能合理模擬超級單體,導致單體過早衰亡,3 km的模擬能更好地模擬超級單體的低層旋轉等特征,1km的模擬能進一步改善對低層旋轉強度的模擬。Sobash等分別使用3 km和1 km分辨率的模式對497個強天氣事件進行預報,評估發現相比3 km分辨率,使用1 km分辨率的預報輸出能診斷獲得更準確的龍卷的相關指標(如低層0~1 km的上升螺旋度),這主要是由于1 km分辨率能更好地預報小尺度對流過程。最近,Squitieri等使用3、1、0.333 km三種分辨率對14個中尺度對流系統的結構進行模擬,其研究也發現1 km分辨率模擬產生的對流冷池堆大小、系統傳播速度、3 h降水率以及9 h降水范圍較為接近實際觀測。水平分辨率對地形降水的影響相對較易理解。在這方面,Heim等通過使用4.4~1.1 km分辨率對阿爾卑斯山區濕對流連續9天的模擬顯示,地形細節對白天的降水分布影響不是很明顯,只是延遲對流的觸發,影響降水量日變化的幅度,但對夜間山谷中對流冷池觸發的降水強度和持續時間影響十分重要。除了確定性預報以外,Schwartz等評估了3 km集合預報和1 km集合預報32 d的降水概率預報結果,發現1 km的集合預報質量通常優于3 km的集合預報(圖2),即雖然二者預報偏差相似,但前者在大多時刻的ETS評分和可測概率(POD)都比后者高,誤報率也比后者低。

圖1 使用分辨率分別為1 km(a)和 0.125 km(b)所模擬的沿颮線深對流區的相當位溫剖面(摘自Bryan等[12]) Fig. 1 Along-line cross sections of equivalent potential temperature from simulations of a squall line using 1 km (a), and 0.125 km (b) grid spacing (from Bryan等[12])

圖2 分別用3 km(黑實線)和1 km(紅實線)分辨率對2013年5月15日—6月15日發生在美國中西部的降水的集合預報試驗結果(摘自Schwartz等[10]) (a)偏差;(b)ETS;(c)可測概率(POD);(d)誤報率(FAR) Fig. 2 Bias (a), ETS (b), POD (c), and FAR(d) of ensemble precipitation forecasts with 3 km (black) and 1 km (red) grid spacing for the central and eastern portions of the United States between 15 May and 15 June 2013 (from Schwartz等[10])

在熱帶氣旋的模擬研究中,對比Liu等對颶風Andrew(1992)的6 km分辨率的模擬,Yau等的研究展現了使用2 km分辨率能更好地模擬颶風Andrew(1992)的螺旋雨帶。同樣,Li等使用1.33 km分辨率能模擬出臺風“龍王”(2005)在登陸時由埃克曼泵吸和非絕熱加熱產生的具有14 h生命史的螺旋雨帶。Gentry等比較了8~1 km對颶風Ivan(2004)的模擬結果,顯示越高分辨率所模擬的颶風最強風半徑越小,眼墻越窄,強度越強。他們還發現格距小于4 km的模擬中開始出現多邊形眼墻(圖3),這可能是因為高分辨率模擬能分辨出眼墻內由切變渦度產生的更小尺度的渦旋結構和演變過程。事實上,一些數值模擬研究顯示高分辨率網格距的使用能使熱帶氣旋的強度達到甚至超過Emanuel所提出的最大理論強度(MPI)。颶風Patrica(2015)的強度創東太平洋熱帶氣旋強度的歷史紀錄,Qin等使用333 m水平分辨率、55垂直層能模擬出其強度和增強速率,模擬的最大地面風達到92 m?s,地面最低中心氣壓達到884 hPa,每小時地面風速增強速率達2.5 m?s,中心氣壓降低速率達4.5 hPa,最強風半徑小于10 km。

圖3 使用分辨率分別為8 km(a),6 km(b),4 km(c),3 km(d),2 km(e)和1 km(f)對臺風所作WRF模擬的 850 hPa位渦(PV),位渦特征與影響臺風強度的對稱或非對稱過程密切相關(摘自Gentry等[22]) Fig. 3 Horizontal distribution of potential vorticity (PV) at 850 hPa from the 8 km (a), 6 km (b), 4 km (c), 3 km (d), 2 km (e), and 1 km (f) simulations of a tropical cyclone. The characteristics of PV are associated with the symmetric and asymmetric processes crucial to TC intensity (from Gentry et al[22])

相較于水平分辨率對數值模擬影響的研究,關于提高垂直分辨率的研究較少。但早在20世紀90年代,已有學者提出水平分辨率和垂直分辨率需要協調一致,模擬預報結果才能更合理。使用與水平分辨率不協調的垂直分辨率,可能會產生更多“噪聲”。目前,很多模擬研究已采用超過70垂直層和幾百帕的模式層頂。Aligo關于垂直分辨率對美國中西部降水預報的研究指出,增加不同高度的垂直分辨率會產生不同的影響,增加低層分辨率有助于邊界層通量模擬,增加融化層以上的垂直分辨率有助提高降水預報。但同時增加不同高度的分辨率可能會由于邊界層過程和云微物理過程共同作用導致降水預報評分并不一定提高。Zhang等對颶風模擬的研究發現增加低層垂直分辨率會增加水汽輻合,影響較高層的潛熱釋放,從而增加模擬颶風的強度。他們認為增加低層垂直分辨率與增加水平分辨率一樣,有利于增強低層輻合。

總之,提高水平、垂直分辨率往往能更好地模擬中小尺度結構特征。但提高分辨率和改善模擬、預報并非簡單的線性關系,還存在不確定的結論,尤其是當分辨率提高到次千米以下。這其中一個重要原因是當前的各種物理方案是依賴于網格分辨率的,不考慮物理方案所包含的假設而簡單地提高模式分辨率是不合理的,在第4節中,我們將針對這個問題進行進一步討論,同時還討論當模式水平分辨率位于積云對流和湍流顯式或隱式處理的所謂“灰色區”時的問題。此外,在不同的地理區域,由于地形復雜程度及其對降水觸發的強迫作用不同(如Heim等),提高分辨率對降水預報的影響也不同。

2 模式初始條件的重要性

模式初始條件的質量對中、小尺度對流系統的可預報性十分重要。實際個例模擬的初始條件通常是從全球分析資料插值得到的。全球分析資料是以全球預報場為初猜場,然后利用多種氣象觀測,包括輻射形式的遙感資料改善初猜場。但是所用到的氣象觀測資料分辨率一般較低,遙感資料則不容易直接用于數值模式,而另一些特別的觀測(如雷達)只能覆蓋小部分區域,存在不同程度的觀測誤差,且動力一致性有所欠缺。一般來說,模擬的最初階段主要依賴于初始場的質量,其中包括觀測和分析誤差。如對在初始或臨近時刻已發生的對流事件,初始場中與之有關的大氣穩定度、濕度、水凝物以及對流尺度環流直接影響到對流系統的大小、強度和未來的發展。如對流系統是在模式起步一段時間后產生,初始條件的質量將決定它產生的時間和地點的可預報性。隨著進一步時間積分,模式動力和物理過程的合理、準確度才顯得越來越重要。另外,模式初始化過程需要將低分辨率的分析場插值到高分辨模式格點上,這存在動力場的“spin up”問題,即產生合理的垂直環流和水汽分布。

資料同化是提高初始場質量,將中小尺度信息融入初始場,縮短“spin up”時間的有效方法,特別是熱帶氣旋的模擬。盛春巖等探討了雷達資料同化和提高模式水平分辨率對一次華北暴雨過程預報的影響,結果顯示18 km分辨率試驗中使用雷達資料同化的預報結果比不使用雷達資料同化的3 km試驗結果要好,但提高分辨率對強降水預報也十分重要。使用經過資料同化改善的較粗分辨率的分析場驅動高分辨率模擬能有效提高高分辨率模擬水平,更好的初始條件也能使提高分辨率至1 km在降水預報方面的優勢更顯著。由此可見,合理運用資料同化能有效改善高分辨率數值模擬。

由于高分辨率數值模擬能分辨包括對流尺度等更多尺度的現象,因此資料同化也必須協調處理不同尺度的特征。對于對流尺度,誤差增長往往會更快,而且還存在非線性等問題。四維變分和集合卡爾曼濾波(EnKF)技術是目前較先進的同化方法,更適合對流尺度同化。

多普勒雷達資料時空分辨率高,能夠反映對流的結構特征,是對流尺度資料同化的重要觀測來源。早期,學者根據雷達反射率估計水凝物含量、潛熱等,然后通過潛熱松弛逼近法(nudging)、云分析等方法進行間接同化。這些方法往往依賴于經驗關系,在對流尺度數值模擬中有一定的局限性。變分方法(圖4和圖5)和EnKF技術較適用于對流尺度同化,采用變分方法或EnKF還能直接同化雷達反射率和徑向風。大量研究采用三維變分(3DVar)直接同化雷達反射率和徑向風。相比3DVar,四維變分(4DVar)能使用到同化時間窗里的觀測,將時間傾向項作為約束,因此可以更好地捕捉對流的演變特征。但4DVar所需計算資源較多。Sun等5應用4DVar技術設計了變分多普勒雷達分析系統(VDRAS),借助該系統可以反演得到與對流相關的三維高分辨率風場、熱力場和微物理場。該系統在北京奧運會期間也有應用。EnKF不需要切線性和伴隨矩陣,它能利用集合預報得到流依賴的背景誤差協方差以及不同變量間的協方差。Dowell等、Tong等展現了使用EnKF方法可以分析得到合理的微物理場、 風場和熱力學變量,在幾次循環同化后風暴可以在模擬中成功再現。他們指出流依賴和動力協調的背景誤差協方差是成功同化的重要因素。Xue等應用EnKF方法,根據雷達徑向風和反射率能成功估計對流風暴的水凝物混合比和數濃度。將由集合預報估計得到的流依賴的背景誤差協方差應用在變分框架中的混合同化方法也有不少進展。

圖4 2002年8月31日0300 UTC臺風茹莎登陸韓國前的雷達反射率(摘自Xiao等[48]) (a)觀測;(b)控制試驗;(c)同化雷達反射率試驗 Fig.4 Radar reflectivity at 0300 UTC 31 August 2002, i.e., prior to the landfall of Typhoon Rusa in Korea (from Xiao等[48]) (a) observation, (b) control experiment, (c) experiment with assimilation of radar reflectivity

圖5 對美國中部2002年暖季6天期間1 h降水量預報所做的6個試驗的FSS評分。CTRL:以Eta模式分析場驅動作為控制試驗;GFS:GFS全球分析場驅動的試驗;3DV_CYC:3小時循環同化常規數據;3DV_RV:同3DV_CYC,并同化雷達徑向風;3DV_RF:同3DV_CYC,并同化雷達反射率;3DV_RD:同3DV_CYC,并同化雷達徑向風和反射率(摘自Sun等[52]) Fig. 5 The 6 experimental-averaged fractions skill score (FSS) of hourly precipitation forecasts over the central United States during a 6-day period in the warm season of 2002 (IHOP_2002). CRTL: initialized by Eta Model analysis; GFS: initialized by GFS analysis; 3DV_CYC: initialized by WRF 3DVAR with 3-hourly cycle without radar; 3DV_RV: as in 3DV_CYC but with radar radial velocity; 3DV_RF: as in 3DV_CYC but with radar reflectivity; 3DV_RD: as in 3DV_CYC but with radar velocity and reflectivity (from Sun等[52])

近年來,其他高分辨率的觀測資料,如地面閃電資料的同化也逐漸受到重視,這些資料能彌補雷達觀測無法覆蓋的區域。如何協調利用好各種新型觀測資料,將高分辨率觀測資料有效同化進模式對提高高分辨率模擬和中、小尺度天氣系統預報的質量十分重要。

3 模式物理過程的作用

強降水由積云對流造成,這涉及到準地轉上升運動、鋒面或地形抬升、邊界層過程、非絕熱加熱/冷卻,輻射和云微物理過程等。這些過程跨越不同的時空尺度。即使使用套網格方法,目前計算機能力也無法包含全部12個量級尺度的大氣運動。所有次網格尺度的物理過程都需通過某種參數化出現在模式中。而不同網格分辨率(即6~8網格距)需使用不同參數化。隨著分辨率的提高,一些原來在較低分辨率模式中不可分辨的現象在高分辨率模式中變得可分辨,需要合理描述各種復雜的物理過程以及它們之間的相互作用。對不同的物理過程,在數值模式中存在一個水平分辨率區間被稱為“灰色區”。在這“灰色區”中,所模擬對象的部分特征是可分辨的(即顯式),但其他特征屬于次網格尺度。完全使用(隱式)參數化方案還是顯式方案就變得模棱兩可(圖6)。 其中比較典型的是積云對流和邊界層及云內的湍流混合處理方案,它們的“灰色區”可粗略認為分別在10~1 km和1000~100 m。其實,本文介紹的大部分模擬工作(除了大渦模擬)所使用的水平分辨率都在這二“灰色區”之中或之一。

圖6 數值模擬中積云對流、大氣邊界層和輻射方案灰色區示意圖(根據 https://www2.mmm.ucar.edu/wrf/users/tutorial/201907/dudhia_physics.pdf改制) Fig. 6 Schematic of “gray zones” for convection, planetary boundary layer, and radiation parameterization schemes (Adapted from https://www2.mmm.ucar.edu/wrf/users/tutorial/201907/dudhia_physics.pdf)

首先是積云參數化方案和微物理方案的使用。水平格距為2~20 km的模擬中,對流不完全是次網格現象,一般來說應同時合理使用積云參數化方案和顯式微物理方案。實際上,目前還沒有發展出適合網格距小于10 km的積云參數化方案。此外,由于不同參數化有不同的閉合假設,使用積云參數化的模式分辨率范圍也有所不同。無論是個例模擬研究,還是對業務模式的評估,關于格距為10 km以下的模擬中使用積云參數化方案對模擬預報會產生正效果還是負效果都無法得出一致的結論。有的模式因使用積云參數化方案能更快觸發對流,同時減輕模擬降水過強的問題,從而使得預報效果更好;但有的使用積云參數化方案后產生過于平滑的結構,預報效果變差。 近年來,學者還嘗試研發自識別尺度(scale-aware)的對流參數化方案,這些方案能從完全的參數化平穩過渡到顯式方案,主要針對水平格距在3~5 km的數值模擬,避免在灰色區完全使用傳統的對流參數化方案或云微物理方案帶來的問題。

高分辨率(2~3 km以下)數值模擬需要使用較完善的云微物理方案。微物理方案可分為兩大類,整體水法(Bulk方法)和分檔法(Bin方法)。Bulk方法假設水凝物粒子的總體分布特征能用一個經驗函數來描述,常用的是Gamma函數,根據函數中可變參數的個數可分為單參方案(僅預報粒子比質量)、雙參方案(預報粒子比質量和數濃度)、三參方案(預報比質量、數濃度、反射率因子或密度)等。典型的雙參微物理方案通常包含水汽、云水、云冰、雨、霰、雪等粒子,描述它們之間的轉化過程,預報各種粒子的混合比以及部分或全部粒子的數濃度(圖7)。Bin方案根據粒子的大小、質量等將粒子分成幾十到幾百檔,建立每類粒子的預報方程和它們之間的轉化關系。再加上使用Bin方案需較高分辨率(1 km以下),使得Bin方案對計算資源要求高。一些研究使用Bin方案研究云-氣溶膠相互作用,以及云微物理-動力相互作用。Bin方案結果也作為對照,用于發展和測試天氣氣候模式中的Bulk方案。但目前更多業務模式和理論研究集中于Bulk方法,探討使用多參微物理方案是否能提高高分辨率模擬的質量。

圖7 典型雙參云微物理方案示意圖(摘自Randall等[80]) Fig. 7 Diagram of a typical two-moment cloud microphysics scheme (from Randall等[80])

不少研究發現,使用雙參方案能夠更好地模擬對流的結構和演變、冷池結構、雙偏振量特征等。例如,Putnam等和Wheatley等闡述了在中尺度對流系統模擬中,由于雙參方案模擬得到較多的雪或冰形成,從對流塔向后輸送到層云區,從而改善了層云降水的模擬。還有一些研究顯示使用三參方案能得到更好的龍卷風暴或雹暴模擬。Dawson等指出使用三參方案能模擬出更接近實況的風暴路徑和冷池特征,而冷池強度對龍卷渦旋的生成和演變十分重要。另一方面,也有研究指出微物理方案中對冰相過程的描述,尤其是關于霰或雹的形成過程的描述,對模擬效果有重要影響。雙參方案對霰或雹的密度和下落速度十分敏感。Adams-Selin等指出微物理方案中關于霰的描述對弓狀回波的模擬有顯著影響(圖8),不同的霰粒子譜分布影響模擬的冷池強度,進而影響氣壓擾動、后向入流和上升氣流的發展等。Tao等和Bae等的研究都發現在微物理方案中加入冰雹對于小雨或者云內過程的模擬有所幫助。當然,上述微物理方案對于不同對流系統或不同分辨率也會有不同效果。如Jensen等使用復雜的冰相微物理方案 (ISHMAEL)以及1.33 km和148 m分辨率模擬湖泊效應產生的對流雪帶。由于后者分辨率高能模擬出較強的上升氣流場,由此產生較多的凇化球狀冰粒和霰粒子,所模擬的雪帶回波強度也更接近于實際雷達觀測。綜上所述,網格距小于1 km的數值模擬應盡可能考慮使用較為完善的云、冰微物理方案。

圖8 使用不同霰粒子截斷參數、密度、下落速度參數的不同(理想)模擬試驗產生弓狀回波時的雷達反射率,粗黑線指示弓狀特征(摘自Adams-Selin等[99]) Fig. 8 Simulated composite reflectivity when bow echo begins to appear in experiments using different intercept values, density and fall speed parameters for graupel. Thick black lines represent bowing segments (from Adams-Selin等[99])

邊界層內動量、熱量和水汽的湍流混合過程對中、小尺度對流系統十分重要。但在網格距1 km以上的模式中往往僅需對湍流垂直輸送參數化。當網格距小于100~200 m時,水平湍流混合開始顯得重要。此外,復雜下墊面、大風或不穩定氣流下的湍渦特征很難通過參數化合理描述。故當模式分辨率達到100~200 m,研究人員使用大渦模擬來直接解析可分辨尺度、較為活躍的高能湍渦,而對次網格尺度湍渦進行參數化閉合處理。但目前常見的幾千米的高分辨率數值模擬無法完全分辨湍渦,因此仍然需要使用邊界層參數化方案。邊界層方案在高分辨率模擬中十分敏感,特別是對熱帶氣旋的模擬,畢竟海面上的潛熱和顯熱是熱帶氣旋發展的主要能源。有研究顯示在可分辨對流的分辨率下使用一些邊界層方案(如MYJ方案)會在邊界層產生冷偏差或濕偏差。

當模式分辨率在幾百米到1 km時,湍渦可以部分被分辨,此時出現對于邊界層而言的灰色區。尤其對于對流邊界層,湍渦尺度可接近1 km,使用常規邊界層方案容易導致誤差。有些學者提出scale-aware邊界層方案,主要思想是根據水平格距,調節由常規邊界層方案計算得到的次網格通量所占總通量的比例。但這些scale-aware邊界層方案大部分針對干對流邊界層,對于濕對流邊界層,涉及到邊界層過程和濕對流的相互作用更為復雜,有待進一步研究。

除邊界層內的湍流交換,對流云內次網格湍流在熱量、水分和其他標量的傳輸中也起著重要作用。在這方面,Sun等開展理想超級單體的大渦模擬試驗,將50 m分辨率的模擬逐步粗格化到4 km,對比這些模擬結果顯示,云內垂直和水平次網格湍流輸送通量在500 m網格就具可比性,并表現出與深對流相似的非局部特征。Kirshbaum通過理想數值試驗來理解在模擬地形對流分別使用對流、湍流灰色區分辨率時的情景,其結果顯示在前者灰色區由于地形產生的水平不均勻流場和加熱場,隱式對流在未達山脊前很遠就發生,時間上也提早激發;而后者灰色區則影響到云內夾帶和對流質量通量的強度。

高分辨率模擬其中一個優勢在于能使用高精度的地形和下墊面信息,這也導致在復雜地形區還需要考慮其他問題,例如在輻射方案中需考慮在云覆蓋率和地形影響下的三維輻射、坡面效應,在處理城市下墊面的冠層模式中需考慮建筑群的分布和高度、道路和非自然加熱等因素。此外,盡管隨著分辨率的提高,可以不再使用原有的參數化方案,但是總會面臨新的更小尺度的次網格過程,如小尺度的三維湍流輸送過程,需要針對這些過程建立新的合理的參數化方案。

4 中、小尺度對流系統的高分辨率數值模擬

由于大多災害天氣與對流尺度環流有關,若要開展災害天氣機理的研究,有必要取得對流系統的高分辨率數值模擬。近年來有不少關于中、小尺度對流系統的高分辨率數值模擬研究。本節介紹使用1 km和次千米分辨率模擬暴雨對流系統、雹暴、龍卷和熱帶氣旋的幾項具有一定代表性的工作, 同時也提及其他相關的高分辨率數值模擬研究。

4.1 暴雨對流系統

梅雨鋒暴雨是影響我國東部的主要災害天氣之一。梅雨鋒暴雨中心的強度和落區受中小尺度對流演變的影響,要深入理解梅雨鋒暴雨的發生發展機制,需要對這些對流及其環境的精細特征進行細致分析,而目前的觀測仍無法提供足夠高分辨率的資料,因此需要依賴高分辨率數值模擬。已有一些研究使用1 km左右,甚至幾百米分辨率的數值模擬結果對梅雨鋒暴雨系統的發生發展機制進行深入分析。Zhang等使用MM5模式,設置四重網格(12/4/1.33/0.444 km),成功模擬了2003年7月4—5日江淮流域中尺度對流系統引發的暴雨過程。四重網格均采用Tao等提出的包含三類冰相粒子的微物理方案和修改的Blackadar邊界層方案,只在最外面兩重使用積云對流方案。該模擬的其中一個關鍵技術在于他們使用外場觀測資料改善NCEP分析場,并且額外使用一個覆蓋范圍遠大于12 km區域的分辨率為36 km的模擬區域,模擬起始時間比12 km區域提前12 h,36 km區域的輸出結果作為12 km區域的初始條件及邊界條件,在36 km區域上使用動力松弛逼近法,以改善提供給12 km區域的邊界條件。

結合觀測資料和高分辨率模擬,他們分析了此次暴雨過程的多尺度特征,包括對流觸發和組織機制。圖9展現了Zhang等中分辨率為444 m的區域模擬的對流系統成熟階段的雷達反射率和局地浮力。其中,計算局地浮力時減去了由32點滑動平均得到的參考虛溫。由于非降水的云水和云冰對雷達反射率的貢獻很小,圖9沒有反映這兩類粒子。如圖所示,分立的(而非連續的)對流單體首先在濕潤的西南季風氣流前觸發,其后在低層輻合區(約10~15 km寬)上是發展高至150 hPa的強降水對流,尾部是變弱的云柱(大部分是層狀降水云)。在550 hPa附近能清楚看到融化層(雷達反射率垂直梯度大值區)。前沿主雨帶由蒸發性下沉冷出流與系統前方高相當位溫氣流匯合而激發生成,其中的上升速度最大可達25 m?s。隨后,這些對流熱塔相對于前沿主雨帶向系統后方移動,并逐漸衰弱,形成層狀云區。在整個對流系統中,包括在層狀云區域,這些不同強度的降水對流熱塔的尺度從2~5 km不等。總的來說,444 m網格的模擬結果更好地重現了對流系統的演變,抓住了多個對流風暴的觸發及其隨后組織形成中尺度對流系統,以及前期對流活動遺留冷空氣造成的中高壓等。他們利用模擬結果進一步闡明了在此次過程中,濕空氣沿著梅雨鋒前的冷空氣頂的等熵抬升對于對流觸發十分重要;夜間西南急流提供充足的不穩定能量和水汽;對流沿著相同路徑形成并移動導致了暴雨的發生。這項研究也展現了使用1 km以下的水平格距可能能更好地模擬對流結構以及相關的β中尺度環流,提高梅雨鋒暴雨的預報水平。

圖9 對2003年7月4—5日發生在江淮流域的一次梅雨鋒暴雨模擬的Vis5D 三維圖像 (a)雷達反射率;(b)局地浮力(紅/藍色陰影分別指示±1℃等值面,綠色陰影表示云,摘自Zhang等[119]) Fig. 9 Vis5D plots of (a) radar reflectivity; (b) local buoyancy from the simulation of torrential rainfall occurring over Yangtze-Huaihe River basin during 4-5 July 2003. Red/blue shadings are for isosurfaces of ±1℃. Green shadings indicate volumes of cloudy air (from Zhang等[119])

需要指出的是上述444 m網格距模擬還使用邊界層參數化,顯然這是屬于邊界層顯式與隱式處理的灰色區(見圖6)。最近,Huang等使用500 m網格距按照WRF大渦模擬方案(即不包含邊界層參數化)模擬了2017年5月6—7日發生在廣州的極端降水事件。他們還與包含邊界層參數化方案的模擬實驗比較,發現后者模擬結果與實際觀測相比不如未包含邊界層參數化的大渦模擬結果。但Yin等使用1.333 km網格和邊界層參數化的WRF也模擬了這一極端降水事件,所模擬的最大降水比Huang等的模擬更接近于實際觀測值。由于這二者模式初始條件、地表處理和其他模式框架都十分相似,這不同的降水模擬結果可能與使用不同的云微物理方案有關。需要指出的是,目前許多使用1 km左右網格距模擬對流觸發、極端降雨以及中尺度對流結構和演變的研究都使用了邊界層參數化;本文前面介紹的大多有關模擬也是如此。不管如何,未來使用1000~200 m網格距來模擬中、小尺度對流系統時有必要開展對邊界層顯式和隱式處理的相對重要性研究。

4.2 冰雹對流系統

冰雹給生命財產安全帶來嚴重威脅,但目前數值模式對冰雹的預報還有一定難度,尤其是地面降雹。此外,相對于其他強災害中尺度對流系統,雹暴的模擬研究工作也較少。Li等針對2005年6月10—11日我國東北地區一次強對流過程,使用WRF模式研究初始環境水汽量對冰雹產生和降水模擬的影響,模擬中設置9/3/1 km三重網格,采用包含霰和雹粒子的雙參Milbrandt-Yau方案,結果顯示在前期,冰雹和降水率隨初始水汽量增加而單調遞增,降水增強造成冷池更強,并維持更長時間,影響深對流發展,使得后期冰雹的產生和初始水汽量呈現非線性變化關系。

Luo等使用雙重(3/1 km)網格的ARPS模式,對2015年4月28日發生在江蘇西南地區的一次多單體雹暴進行模擬。他們對比了使用單參、雙參和三參數版本的Milbrandt-Yau方案在模擬冰雹最大可能直徑、尺度分布、累積質量、數濃度等方面的差異。與S波段多普勒雷達觀測結果比較顯示,使用三參方案的模擬效果較好(圖10)。不同方案對冰雹粒子譜分布的形狀參數描述的差異導致雹生長過程明顯不同,加上粒子分選的差異,共同導致模擬結果的差異。單參方案沒有粒子分選過程;在雙參方案中使用固定的形狀參數(設為0)會導致過度的粒子分選,使得冰雹的粒子譜分布偏向于大雹;而三參方案能獲得更合理的形狀參數,產生更合理的粒子譜分布,模擬的最大可能雹直徑、地表累積雹濃度等優于單參和雙參方案。由此可見,在高分辨率云模式中使用合理的形狀參數對于冰雹各微物理參數的模擬和冰雹事件的預報十分重要。

圖10 2015年4月28日發生在江蘇儀征到常州、無錫和蘇州一帶的多單體雹暴觀測(a,OBS)和用ARPS 模式中三參微物理方案(b,CNTL)預報試驗結果所診斷的最大雹直徑(陰影,mm)分布(摘自Luo等[96]) Fig. 10 Maximum estimated sizes of hail from observation (a, OBS), and simulation with the three-moment microphysics scheme in ARPS model (b, CNTL) for a multicellular convective system occurring over Yizheng, Wuxi, and Suzhou on 28 April 2015 (from Luo等[96])

后來,Labriora等使用500 m網格距評估了ARPS模式中Milbrandt-Yau雙參和三參數方案以及NSSL多密度冰融化雙參方案對美國一次冰雹事件的0~90 min預報,這三方案都能報出大致的冰雹范圍,但預報強冰雹發生范圍的能力較低。最近Mansell等在NSSL三參數方案中加入允許改變分布形狀并適應融化對冰雹和雨滴尺寸分布的一些直接影響,能模擬出冰雹粒子融化時由于重力作用使得冰雹和雨滴大小重分布的狀況。

4.3 龍卷結構

由于龍卷的尺度小、維持時間短,我們對龍卷的觀測研究仍十分困難。目前對龍卷渦旋結構的理解大多來自數值模擬結果的分析。自20世紀70年代末Klemp等使用2 km格距和(

x

,

y

,

z

)維度為24×24×20的格點模擬龍卷型超級單體以來,到近10年來隨著超大型并行計算機性能的快速提升,人們已能使用20~100 m分辨率來模擬龍卷系統的結構。雖然大多數研究仍是基于探空廓線的理想模擬,模擬區域也較小,但所取得的模擬結果可以提供目前觀測資料難以解析的中、低層三維渦旋風場和熱力結構。如Orf等使用Cloud Model 1(CM1)變網格以及雙參方案模擬了2011年5月24日在美國俄克拉何馬產生EF5級別龍卷的超級單體。在模式區域中心使用了水平和垂直格距為30 m,(

x

,

y

,

z

)維度為2200×2200×380,即近20億的格點數; 時間步長為0.2 s。變網格模式中粗網格三維區域范圍為160 km×160 km×20 km。圖11展現了龍卷生成時的三維流向渦旋。所示的流向旋流(SVC)和Schenkman等模擬的水平旋流相似,但SVC能延伸到前側更深處,而且與水平旋流不同,由于自由滑動的下邊界條件,SVC渦度不能由表面摩擦產生。在SVC傾斜至上升氣流前,SVC的水平渦度由前側下沉邊界(FFDB)相關的浮力梯度產生。隨著FFDB的浮力梯度增大,SVC增強。在龍卷發展前的10~12 min,SVC顯著增強。隨著SVC增強,SVC傾斜同時加強了低層中尺度渦旋的旋轉以及隨后的上升氣流。上升氣流的加強和旋轉引起的向上的擾動氣壓梯度加速一致。加強的上升氣流導致其下方形成輻合區。上述龍卷結構能在模式中傳播120 km、維持2.5 h以上。

圖11 CM1模式在5100 s時模擬的相對于龍卷的三維流向渦旋:從左側垂直風切變引起的低層水平渦旋到龍卷中心的渦旋塔(摘自Orf等[139]) Fig.11 Storm-relative 3D streamwise vorticity (from the low-level horizontal vorticity induced by vertical wind shear on the left to the vorticity tower in the center of tornado) at t = 5100s of the simulation by CM1 (from Orf等[139])

Sun等使用WRF模式,設置1~5層嵌套網格,最內層網格分別為4000、1333、444、148、49 m,嘗試模擬2016年6月23日江蘇阜寧EF4級龍卷。結果顯示,使用444、148和49 m分辨率能分別模擬出強度達EF1、EF2和EF3的龍卷。49 m分辨率的模擬結果(圖12)中,多個“抽吸渦旋”沿著高渦度環發展,從而形成具有多渦旋結構的龍卷。次渦旋、主渦旋環流和系統移動速度的疊加產生局地強風區。根據他們的研究結果,在業務預警應用中需要使用至少500 m分辨率才能抓住龍卷渦旋特征,而研究龍卷動力結構,則需要采用至少50 m分辨率的模擬。

圖12 使用WRF模式49m分辨率模擬江蘇阜寧EF4級龍卷的結果:(a)地面(z=10 m)水平風(填色,m?s–1)和垂直渦度(粗線,s–1);(b)最大地面風速軌跡帶(摘自Sun等[142]) Fig. 12 Results from the WRF simulation of the Funing/Jiangsu EF4 tornado with horizontal grid spacing of 49m: (a) surface (z =10 m AGL) horizontal wind vectors and speed (shaded, m?s-1), and vertical vorticity (thick contours, s-1); and (b) swaths of the maximum surface wind speed (from Sun等[142])

4.4 熱帶氣旋的精細結構

東亞和北美洲的不少國家受熱帶氣旋影響嚴重。雖然熱帶氣旋是α中尺度的系統,但其內部復雜的結構特征可能引發局地性的暴雨、大風等災害。為更好地滿足防災減災的要求,需要深入研究熱帶氣旋的精細化結構。在這方面,已有不少使用1 km和次千米網格距模擬熱帶氣旋生成、爆發性增長、極端降水和精細結構的工作。由于篇幅有限,在此僅介紹兩篇揭示不同臺風精細結構的大渦模擬研究。

Ito等利用日本氣象廳區域天氣預報模式的大渦模擬版本對臺風整體進行大渦模擬(圖13)。他們的模擬沒有使用嵌套網格,避免了不同分辨率網格之間的邊界問題。模式分辨率設為100 m,可分辨臺風邊界層的大渦,沒有使用邊界層參數化方案;設置60個垂直層,高度為23 km,其中538 m以下的垂直層間隔小于100 m。100 m網格的初始場是由使用同樣區域,但分辨率為2 km的網格模擬輸出插值得到的。該模擬使用了9216個計算節點,占了所用的K超算系統將近1/8的資源,積分1個小時需要約9.5 h,臨時輸出文件大小達17 TB。分析模擬結果顯示,大渦模擬的近地面最大風速和2km分辨率網格模擬的相近,但是大渦模擬呈現出更精細的邊界層結構特征,主要涉及三種類型的卷流。A型卷流主要位于大風半徑外;B型卷流靠近大風半徑;C型卷流位于大風半徑內,可導致逆梯度輸送動量,引起強陣風。該模擬研究應是首個發現B、C型卷流的研究。雖然圖13展示的是作者們使用理想初始條件所作的數值試驗,他們也使用類似的模式框架對2019年9月給日本造成嚴重損失的第15號臺風法茜進行了10 h整體大渦模擬,在美國氣象學會100周年年會期間作了報告(https://ams.confex.com/ams/2020Annual/videogateway.cgi/id/517686?recordingid=517686),其動態影像十分吸引眼球,引起人們強烈的科學興趣。

圖13 (a)臺風整體大渦模擬的云的三維視圖;(b)圖13a臺風中心附近的放大圖,橙色方框指示C型卷流的區域(摘自Ito等[153]) Fig. 13 (a) Three-dimensional view of simulated cloud amount from the large eddy simulation of entire tropical cyclone; and (b) is a close-up view of the typhoon inner core region in (a). The orange rectangle indicates the region of Type-C rolls (from Ito等[153])

大渦模擬可以為深入研究熱帶氣旋的精細化結構提供極高分辨率的信息。以往觀測發現熱帶氣旋邊界層存在幾百米到幾千米尺度的小尺度系統,甚至龍卷尺度渦旋,這些渦旋的垂直渦度峰值和水平尺度與弱龍卷相當,有時伴隨著10~30 m?s的極強上升運動,并可能造成局地強陣風。目前對熱帶氣旋邊界層的觀測還較少,精度也不足,難以對龍卷尺度渦旋開展精細分析,需要借助百米級的超高分辨率模擬。Wu等使用 6重(9/3/1/0.333/0.111/0.037 km)嵌套網格和75個垂直層(其中2 km以下有19層)的WRF模式模擬研究臺風邊界層的精細結構。在分辨率為1 km以下的網格使用大渦模擬技術,其余網格使用YSU邊界層參數化方案。他們成功模擬出和觀測較一致的龍卷尺度渦旋特征(圖14),包括強上升速度和近地面局地大風等,龍卷尺度渦旋附近風速可達70 m?s,并發現龍卷尺度渦旋普遍存在于臺風眼墻的內測,它們很可能影響近地面極端陣風的分布。

圖14 模擬龍卷尺度渦旋相關的擾動風場。流線是水平擾動風,暖(冷)色流線表示上升(下沉)運動;箭頭指示近地面風,填色是風速。臺風中心在左上角方位(摘自Wu等 [159]) Fig. 14 The streamlines of the horizontal perturbation winds associated with the simulated tornado-scale vortex. The warm (cold) color of the streamline denotes the upward (downward) vertical motion. Vectors represent the near-surface winds with the wind speed shaded. The typhoon center is located far outside from the left corner (from Wu等 [159])

5 未來挑戰

綜上所述,目前中、小尺度對流系統的數值模擬研究開始進入次千米分辨率的階段。提高模式分辨率有助于減小由物理過程參數化造成的不確定性,從而改善對中、小尺度對流系統強度、結構和演變的模擬,整體提高數值天氣預報水平。更高分辨率的模擬結果也有助于加深對這些對流系統及其相關天氣現象的形成機制的理解。在不久的將來,越來越多的業務天氣預報區域模式會使用1 km左右的分辨率,我們也會看到更多中、小尺度對流系統真實個例的次千米級模擬,甚至大渦模擬。此外,取決于計算能力的提升程度,1~4 km對流分辨率的全球預報模式也可能在不久的將來投入業務運行。我們相信,到那時通過多重嵌套或變網格技術將會逐漸淘汰區域模式,從而不斷增加使用高分辨率、不存在側邊界問題的一體化全球模式。

但是,未來提高模式分辨率以改善對中、小尺度對流系統的數值預報、模擬仍會面臨不少挑戰。首先,如前所述,目前高分辨率模擬都是通過多重嵌套網格實現的,很多情況下僅有部分對流系統被細網格覆蓋。增加細網格覆蓋面積或提高模式分辨率,計算量將呈幾何級數增長。本文中介紹的幾篇大渦模擬工作幾乎已將能獲得的計算能力推至極限,所消費的經費也是十分“昂貴”。由于增加時鐘頻率需要耗費昂貴的電力,目前而言時鐘頻率已難以大幅增加。從2000年左右起,超級計算機大多通過增加并行核數來提高計算能力。處理器之間的信息交換速度相對較慢,對數值預報模式中所發展的任何新程序并行化,以及重構模式代碼使能高效利用大量的處理器,這些問題需要未來計算方式和能力有革命化的發展(如未來量子計算機的推廣運用)。

其次,隨著觀測數據(包括遙感數據)越來越多樣化,數據量越來越大,需要發展更先進、有效的客觀分析和數據同化方法以充分利用觀測資料,獲得更高質量的動力一致的模式初始場。使用融合了豐富觀測信息的高分辨率分析場來驅動高分辨率數值模擬對改善模擬效果有重要影響,尤其是提高小尺度對流系統(如龍卷風暴,由不均勻地表條件所引發的局地風暴等)初始形成的時間和地點的可預報性。但是,我們也需要明白,由于無法準確得知模式的最優初始狀態,模式物理方案也存在很多不確定性,我們難以期望獲得對中、小尺度對流系統的完美模擬。

由于模式物理方案具有尺度依賴性,自識別尺度的概念近年來也得到越來越多的關注。在理想情況下,自識別尺度的模式物理方案應該能夠精確地從滿足靜力平衡的隱式方案(即大于10 km)為主,自然過渡到模擬(非靜力)對流云的顯式方案(即1~2 km),再到大渦模擬(即10~100 m)的網格距尺度。由于我們對不同尺度物理過程的認識還十分有限,要發展如此理想的自識別尺度物理方案還是十分困難的。另一方面,即使在較高分辨率的模式中也存在難以分辨的現象,也需要發展有關的合理的參數化方案。如在次千米分辨率中,一些對流要素仍無法被分辨;在幾十米分辨率的情況下,一些云微物理過程也還需要進行參數化描述。獲得合適的觀測以驗證這些對流云過程也仍是十分艱巨的任務。在模式中描述邊界層內和積云對流中的三維湍流混合輸送過程(特別是大風條件下)也面臨著類似的挑戰。

本文提及許多高分辨率數值模擬展現了目前觀測手段還無法得到的對流系統內核結構。這些模擬的小尺度對流結構仍需由高分辨率的觀測驗證,才能真正確信其模擬結果。而形成這些結構的有些物理過程(如潛熱釋放)是目前無法測量的或可獲得的相關氣象信息甚少,以至于難以提供理論上的解釋。顯然,未來需要發展合適的觀測手段,獲得更精細的觀測資料來幫助驗證模擬結果的真實性。

正如大多運行模式者所知,每一個中、小尺度對流系統個例的高分辨率數值模擬都將產生動力學一致的海量資料(如Orf等,Ito等和Wu等的大渦模擬)。但以模式模擬為主已發表的很多文章僅利用了其中很有限或僅對作者有用的信息,很多寶貴的有科學意義的氣象信息(特別是一些讀者感興趣的資料)就被“浪費了”。如何共享這些模式資料并充分診斷分析、理解其所提供的中、小尺度信息十分重要。這需要發展更有效的多維診斷分析軟件(如本文引用的幾張三維圖)。通過“人工智能”方法來提高對模式資料診斷分析的能力可能是一個好途徑。此外,還需加強科技人員之間的合作,以最終將我們對模擬結果的充分理解上升為中、小尺度對流系統可預報性、中尺度動力學和降水物理學的理論。

總之,高分辨率數值模擬和預報已得到廣泛應用,但仍需從初始條件、物理方案等多方面改善模擬預報的效果,尤其是當分辨率進一步提高到1 km乃至幾百米時,還有不少問題有待深究。例如,需要配合發展合適的對流尺度同化方法以提供與模式分辨率相協調的更優的初始條件,需要改進、合理使用顯式微物理方案和邊界層模擬技術等。與此同時,未來需要更多精細的觀測來進一步驗證高分辨率數值模擬結果,也為改進物理方案提供更多信息,使物理方案適用于更高分辨率的模擬。

主站蜘蛛池模板: 国产亚洲视频免费播放| 99re免费视频| 老司机午夜精品视频你懂的| 人人91人人澡人人妻人人爽| 国产成人三级| 欧美另类第一页| 一区二区三区在线不卡免费| 欧美日本不卡| 国产日韩欧美精品区性色| 视频国产精品丝袜第一页| 熟女视频91| 国产自产视频一区二区三区| 亚洲国产中文精品va在线播放| 国产女人水多毛片18| 国产喷水视频| 青草视频网站在线观看| 在线观看亚洲成人| 伊人久久综在合线亚洲2019| 欧美日韩国产高清一区二区三区| 精品一区二区三区水蜜桃| 亚洲中文字幕久久精品无码一区| 免费无码AV片在线观看中文| 国国产a国产片免费麻豆| 欧美日韩激情| 丁香婷婷久久| 好紧太爽了视频免费无码| 中文字幕亚洲综久久2021| 午夜性刺激在线观看免费| 91精品专区| 日韩av无码DVD| 91九色视频网| 精品国产自在在线在线观看| 一本大道AV人久久综合| 欧美午夜精品| 国产h视频在线观看视频| 热九九精品| 在线亚洲精品自拍| 久久一级电影| 91成人试看福利体验区| 亚洲综合片| 伊人色天堂| 国产精女同一区二区三区久| 欧美区日韩区| 精品视频一区二区观看| 1769国产精品视频免费观看| 亚洲精品欧美日本中文字幕| 成人午夜久久| 少妇精品在线| 亚洲日韩第九十九页| 日韩欧美中文字幕在线韩免费 | 亚洲日本精品一区二区| 中国成人在线视频| 五月丁香在线视频| 亚洲综合一区国产精品| 国产成人精品在线1区| 日韩精品无码一级毛片免费| 色婷婷亚洲综合五月| 亚洲精品在线影院| 成人免费一区二区三区| 嫩草在线视频| 爱色欧美亚洲综合图区| 久久精品中文无码资源站| 538国产在线| 久久精品视频一| av免费在线观看美女叉开腿| 久久不卡精品| a毛片免费看| 精品国产自在现线看久久| 国产麻豆福利av在线播放| 亚欧成人无码AV在线播放| 亚洲AV无码乱码在线观看代蜜桃 | 人妻无码中文字幕一区二区三区| 台湾AV国片精品女同性| 亚洲永久视频| 无码AV日韩一二三区| www.亚洲天堂| 伊人久久影视| 天天躁夜夜躁狠狠躁躁88| 麻豆精品久久久久久久99蜜桃| 国产成人综合亚洲欧美在| 在线无码私拍| 国产aⅴ无码专区亚洲av综合网|