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

基于DART +WACCM模式搭建的平流層同化、天氣預報和氣候預測模型研究

2022-12-03 02:40:30謝飛田文壽鄭飛張健愷陸進鵬
大氣科學 2022年6期
關鍵詞:大氣

謝飛 田文壽 鄭飛 張健愷 陸進鵬

1 北京師范大學系統科學學院, 北京 100875

2 蘭州大學半干旱氣候變化教育部重點實驗室, 蘭州 730000

3 中國科學院大氣物理研究所國際氣候與環境科學中心, 北京 100029

4 南京大學大氣科學學院,南京 210023

1 引言

平流層位于大氣10~50 km之間的高度層,屬于臨近空間的中低層部分。平流層空氣的垂直運動非常弱,大氣演變過程比較緩慢,穩定度高,也幾乎沒有對流層中的那些復雜且瞬變的天氣現象。Baldwin and Dunkerton(2001)通過對日平均NCEP(National Centers for Environmental Prediction)再分析資料的分析證實,平流層中10 hPa以上的強位勢高度異常會隨著時間向下傳輸,并在15天左右傳輸到對流層中,影響對流層的天氣氣候變化。近些年來,國內一些學者的研究甚至發現,我國所在的東亞地區,其氣候變化也會受平流層極渦、熱帶平流層過程的影響(例如:胡永云, 2006; Chen et al., 2011; 田 文 壽 等, 2011; Ren and Hu, 2014)。一方面,由于平流層環流異常可以對對流層天氣過程與氣候變化造成不可忽視的影響;另一方面,平流層是未來臨近空間低動態飛行器活動的新區域。研究平流層的環流、溫度及與其相關的化學成分變化,既能進一步理解對流層氣候變化,也能為臨近空間低動態飛行器飛行活動提供氣象安全保障。因此,平流層大氣環流、溫度變化及與其相關的化學成分變化已受到了國內外學者的廣泛關注(例如:Yi et al., 2002; 呂達仁和陳洪濱, 2003; Zheng et al.,2004; 陳 洪 濱 等, 2006; Li et al., 2007; Chen et al.,2008; 陳月娟等, 2009; 劉毅和劉傳熙, 2009; 卞建春等, 2011; Xu et al., 2015)。

平流層作為傳統大氣科學和空間科學研究的過渡地帶,前人主要集中在平流層動力—物理—化學過程機制方面的研究,而對平流層環流、溫度和化學成分變化的預報和短期氣候預測方面的相關工作卻非常少。在平流層大氣的研究中,一方面,平流層觀測數據少、種類少、缺失值多,而大氣環流模式對平流層環流隨時間變化的模擬水平與觀測資料之間還有一定差距。數值同化技術可能為平流層研究產生更接近觀測值的、更完整的數據集。另一方面,由于平流層與對流層的相互作用,數值模式中平流層變化的預報預測效果可能會對對流層整個動力過程和熱力場的預報預測,甚至對全球天氣預報和短期氣候預測結果都有極其重要的影響。天氣氣候變化的預報預測技術的高低,受數值模式本身包含的化學反應和物理過程完整度、物理參數化方案優化程度、邊界條件等影響;另外,預報預測數值試驗的初始場精度也是影響預報預測結果的一個重要條件。當下,使用觀測資料進行數值同化就是估計數值模式初值這一問題的重要方法和手段。21世紀以來,平流層大氣探測和反演技術在不斷的發展(呂達仁和陳洪濱, 2003),以往同化對流層氣象觀測資料的技術正在向平流層大氣延伸。

針對平流層大氣資料數值同化技術,近些年國內外的一些學者對此核心問題開展了研究。比如,國際上就有學者基于三維變分同化方法開發了中高層大氣數值同化技術。Jackson(2007)通過英國氣象局資料同化系統(基于三維變分同化方案)將MLS(Microwave Limb Sounder)衛星的臭氧觀測資料同化到了數值模式中。他們的研究中發現,對MLS觀測資料同化后,模擬結果中平流層臭氧體積混合比的誤差標準差和平均分析誤差降低了。Pierce et al.(2007)指出將太陽掩星臭氧體積濃度觀測資料(Stratospheric Aerosol and Gas Experiment III)同化到模式中,能夠顯著提高RAQMS(Realtime Air Quality Modeling System)模式中低平流層和上對流層的臭氧體積濃度的模擬效果。肖存英等(2017)也使用三維變分同化方案,在模式中連續同化了全球10~96 km(臨近空間)高度范圍內MLS衛星臭氧數據。此外,三維變分同化方法還被應用到了準兩日振蕩(McCormack et al., 2009)、極地云(Siskind et al., 2011)和平流層爆發性增溫事件(Coy et al., 2009; Ren et al., 2011; Wang et al.,2011)的研究中。四維變分資料同化,則是在三維變分資料同化的基礎上,加入時間維度信息,把多個時刻的大氣狀態數據聯系起來,找到數據之間的變化關系。通過四維變分方案,Errera et al.(2008)在一個化學傳輸模式中,同化了臭氧體積濃度觀測資料。他們的研究發現同化臭氧濃度觀測資料能夠較好地提高模式中平流層臭氧洞結構變化的模擬效果。后來,Elbern et al.(2010)通過結合流依賴模型與四維變分同化技術,研發出了SACADA(The Synoptic Analysis of Chemical Constituents by Advanced Data Assimilation)系統。這是一種新型的中高層大氣化學數值同化系統,此系統甚至可以同時同化觀測的O3、N2O、CH4、HNO3、NO2和水汽等資料。然而,使用三維變分或者四維變分方法進行數據同化也并非沒有潛在的問題。特別是,由于錯誤的設定背景誤差協方差,可能會導致出現虛假的相關關系(Polavarapu et al., 2005)。

集合數據同化是可用于大氣變化研究的另一種數據同化方法。集合數據同化會從集合成員中直接獲得背景誤差協方差,這樣就不需要指定背景誤差協方差,從而降低了由于指定背景誤差協方差而產生的虛假相關性的可能性。因此,集合數據同化方法在大氣模型數據同化中具有一定的潛在優勢。由美國大氣研究中心(NCAR)開發的DART(Data Assimilation Research TestBed)同化工具,采用的集合數據同化方法就是集合卡曼濾波方法。Pedatella et al.(2013)利用DART研究了WACCM(Whole Atmosphere Community Climate Model)模式在同化觀測資料以后,對中高層大氣中的日變率模擬性能的改善。他們發現,在同化了來自輻射傳感器和飛機的溫度和風、衛星的風、低層大氣的COSMIC(Constellation Observing System for Meteorology,Ionosphere, and Climate)折射率,以及中高層大氣的SABER(Sounding of the Atmosphere using Broadband Emission Radiometry)溫度觀測資料后,模擬的全球大氣環流與觀測資料之間的均方根誤差減少了40%。他們指出,WACCM+DART可以重現中高層大氣的日變化規律,能夠深入研究從地表到平流層的真實日變化特征。隨后,Pedatella et al.(2013)利用WACCM+DART重現了2009年平流層爆發性增溫的時空演變過程。他們的研究中,同化了中高層大氣的MLS和SABER衛星溫度觀測資料、COSMIC掩星觀測資料和低層大氣的飛機、探空、衛星云軌跡風等觀測資料。對2009年1月平流層爆發性增溫的模擬顯示,在升溫過程中,平流層的行星波活動顯著增加,這可能是由于行星1波和2波之間的非線性波流相互作用所致。他們通過對全大氣數據同化系統的模擬分析還發現,對流層和平流層之間的大尺度動態耦合過程,對熱層和電離層的變化具有提前數天的潛在預報能力。Sassi et al.(2013)同樣使用WACCM加數據同化技術,研究了2009年1月與2月之間發生的平流層大面積升溫事件,并成功地重現了事件發生后異常環流持續數周的現象,并重點研究了200公里(低層熱層)以下大氣環流的動態響應。

在前人的研究中,他們主要關心的是中高層大氣同化技術的發展,平流層天氣尺度預報和短期氣候預測相關技術的研究涉及較少。而且他們使用的中高層大氣模式和數值同化工具版本相對較老,有必要嘗試使用最新版本的模式和同化工具搭建的模型重新研究這些問題。國內平流層大氣的數值同化研究起步較晚。基于WACCM+DART,敬文琪等(2019)將2016年2月一次平流層爆發性增溫(SSW)過程作為個例模擬進行了臭氧觀測同化試驗,發現WACCM+DART臭氧分析場能夠較真實反映SSW期間北極上空平流層臭氧廓線隨時間的演變特征。但是,研究過程中,只同化了臭氧變量。特別是對于平流層大氣預報和預測方面的技術還沒有受到國內科學家的足夠重視。為了將我國在平流層同化、預報和預測領域的研究推進到國際領先地位,建立基于多源觀測的平流層大氣環境狀況數據庫,開發面向平流層大氣高影響天氣的專業模式和數值預報,既能加強對平流層大氣的基礎認識,推動平流層大氣探測技術的發展,重新認識平流層大氣中關鍵天氣、氣候要素和化學成分在不同時間尺度上的空間分布特征以及平流層大氣與對流層的動力耦合機制,又能解決國家國防活動重大需求,為臨近空間大氣飛行器提供氣象環境安全保障。這個研究中,我們基于WACCM最新的版本WACCM6大氣化學氣候模式和DART最新版本DART_Manhattan同化工具搭建了一個具有多變量數值同化、0~30天天氣尺度預報和30~60天短期氣候預測能力的模型。為了檢測模型的性能,研究中進行了同化、預報和預測3個方面的檢測。以2020年3、4月北極平流層出現的大規模損耗(Hu,2020; Witze, 2020)為主要研究個例,檢測了同化中高層衛星資料后,模型對臭氧洞事件的重現能力。以2020月4月底的同化結果分析場為初值,對5、6月的平流層大氣狀態進行了回報試驗,并對比了預報預測結果與其他再分析資料的差異。文章的結構如下:第2章詳細介紹了同化、預報和預測模型,第3章利用本模型重現了2020年3、4月北極平流層臭氧大規模損耗事件,第4章分析了模型對北極平流層環流和溫度的同化和預報、預測能力,第5章診斷了模型對全球平流層溫度和環流模擬、預報和預測能力的改善,第6章為討論和結論。

2 數據、模式與試驗設計

2.1 數據

MERRA2 的溫度場和緯向風再分析資料,用于同化效果的檢驗。MERRA2是MERRA的第2版。其水平分辨率為0.625°×0.5°。MERRA2的垂直分辨率在中上層平流層約為2~4 km,在對流層上層和下層平流層約為1~2 km,從1000~0.1 hPa共72層。

MLS的臭氧體積濃度衛星觀測資料,用于同化效果的檢驗。MLS是搭載在由NASA于2004年7月15日發射升空的Aura 衛星上的微波臨邊探測器。MLS探測器臨邊探測廓線水平間隔約為200~550 km,垂直分辨率為2.5~5.5 km,空間覆蓋范圍為82°S~82°N。本研究采用的MLS v4.2x標準臭氧體積濃度產品是由240 GHz的輻射值與115 GHz 輻射值獲取的大氣溫度產品反演得到的。根據MLS v4.2x 2 級產品質量和描述文檔(Livesey,2015),261~0.02 hPa范圍之外的臭氧體積濃度資料不推薦用于科學研究用途,因此本研究只使用氣壓范圍為261~0.02 hPa 的臭氧體積濃度資料。

TIMED/SABER臭氧、水汽體積濃度和溫度衛星觀測資料,用于中高層大氣數值同化,SABER是美國航天局TIMED衛星上的四個儀器之一。SABER儀器的主要目標是為我們研究平流層、中間層和熱層低層的能量、化學、動力學和傳輸的基本過程和特征提供必要的數據。SABER通過使用10個通道寬帶的臨邊掃描紅外輻射計對大氣層進行全球掃描,其光譜范圍為1.27 μm至17 μm。SABER可以提供動力學溫度、壓強、位勢高度、痕量氣體O3、CO2、H2O、[O]和[H]的體積混合比等變量的垂直剖面圖。

2.2 WACCM模式

WACCM是美國大氣研究中心近年來發展的“對流層—平流層—中間層—熱層底部”一體化的全球三維大氣模式,采用有限體積動力框架(Lin,2004)。WACCM中的化學過程通過三維全球大氣化學輸送模式MOZART來實現(Kinnison et al.,2007),包含50余種主要的化學成分以及多種中間層與熱層下層區域中的主要離子成分,并考慮了多種痕量氣體的化學反應過程(Garcia et al., 2007)。

本研究采用的模式版本為WACCM的最新第6代版本(WACCM6),它代表了中層大氣直至中間層模擬的最新水平。WACCM6的水平分辨率使用0.9°(緯度)×1.25°(經度)。在垂直方向上,WACCM6 從邊界層向上擴展至熱層底部(5.96×10?6hPa;140 km)共70個垂直層次,垂直分辨率在對流層約為1.1 km,平流層低層1.1~1.4 km 之間,平流層頂(約50 km)高度附近為1.75 km,65 km 以上垂直分辨率約為3.5 km。與以往的WACCM版本相比,WACCM6再現的中層大氣中的溫度、風和微量成分(如水汽和臭氧)的氣候態特征與觀測更接近,能夠更準確模擬SSW導致的平流層變化,改善了全球溫度和平流層環流對火山噴發的響應,對QBO和中層大氣的長期趨勢有很好的模擬水平。其模擬的溫度、風和水汽與觀測相比偏差很小,比以前版本的WACCM小得多。盡管模擬的南半球極地平流層在春季和初夏與觀測相比存在溫度偏差,但WACCM6能夠重現20世紀和21世紀臭氧的演變特征。另外,熱帶地區的模擬和觀測相比還存在一些偏差,這可能是由于熱帶上升流的速度模擬不好造成的。WACCM6對平流層模擬的改善,改善了全球氣候變化的模擬,甚至改善了地表氣候的模擬。WACCM6模擬的高緯變率,特別是冬季海平面氣壓的標準差,與觀測結果更為一致,而且在WACCM6中模擬的阻塞頻率也與觀測更接近。

2.3 DART 同化系統

DART是由NCAR資料同化研究部門研發的一款基于集合EAKF(Ensemble Adjustment Kalman Filter)濾波理論的、適用于資料同化教學、研究和開發的開源軟件工具(Anderson, 2010)。最小二乘框架下的EAKF基本算法(Anderson, 2003)分為兩步:第一步,分別在單個觀測位置上利用標量集合濾波方法,計算每個觀測變量先驗估計集合成員的更新增量;第二步,利用觀測變量先驗估計集合的更新增量和觀測變量先驗估計集合與每個模式狀態變量集合的線性回歸關系,線性化地計算每個模式狀態變量集合的更新增量。集合EAKF相對傳統的EnKF方法的主要優點在于:一方面它不需要通過矩陣運算構建線性調整矩陣,從而耗費大量計算資源;另一方面它將觀測向量分為多個單一觀測標量,分別計算觀測增量,便于并行化,更符合實際的工作需求。目前,DART工具很多觀測類型的同化接口有待進一步開發。

本研究使用的DART同化工具版本為最新的Manhattan版本。DART采用模塊化編程方法,應用集成卡爾曼濾波器,該濾波器的作用是將模型值調整為與一組觀測數據信息更一致的狀態,這需要運行DART輸入的多個觀測實例來生成一個狀態集合。也就是將適合被同化的觀測算子類型應用于每個狀態,以生成模型對觀測的估計,然后將這些估計及其不確定性與觀測值及其不確定性進行比較,最終對模型的狀態進行調整。

2.4 同化、預報和預測試驗設計

本研究一共設計了四組試驗。第一組為歷史重現試驗,模擬時期為2020年3月1日至4月30日。采用WACCM6+DART的同化模擬,同化了SABER衛星溫度、臭氧和水汽資料,溫室氣體強迫 來 自CMIP6(Coupled Model Intercomparison Project Phase 6)2020年3~4月多模式集合平均結果,海溫驅動資料來自Hadley觀測資料。第二組也為歷史重現試驗,模擬時期同樣為2020年3月1日至4月30日。采用WACCM6單獨模擬,溫室氣體強迫來自CMIP6 2020年3~4月多模式集合平均結果,海溫驅動資料來自Hadley觀測資料。其與第一組試驗相比,只是沒有進行數值同化處理。第三組為預報預測試驗,采用WACCM6模型,回報時段為2020年5月1日至6月30日。溫室氣體強迫來自CMIP6 2020年5~6月多模式集合平均預測結果,未來海溫海冰強迫資料來自NCEP CFSv2(Climate Forecast System Model version 2),預報初值為第一組試驗輸出的4月30日分析場。第四組試驗與第三組完全類似,只是試驗初值為第二組試驗輸出的4月30日分析場。四組試驗的模式時間積分步長設置為1800秒。同化試驗中,同化周期為6 h,同化窗區為±1.5 h。四組試驗中,每組試驗5個集合成員,文章中的結果均為5個集合試驗平均結果。表1更清晰地顯示了四組試驗的設計方案。

表1 四組試驗設計Table 1 Design of experiments

第一組試驗在本研究中有2個作用:(1)與第二組沒有進行同化的試驗對比,以說明同化衛星資料以后對歷史重現試驗的改善效果;(2)第一組試驗的2020年4月30日的分析場為第三組WACCM6模式預報預測試驗提供初始場。而第二組試驗的2個作用是:(1)與第一組試驗結果對比(同上);(2)第二組試驗的2020年4月30日的分析場為第四組WACCM6模式預報預測試驗提供初始場。第三組和第四組為歷史回報試驗。第三組試驗的作用是檢測利用WACCM6+DART同化模型為WACCM6提供了預報初值后,WACCM6預報預測模型對歷史過程的回報效果。第四組試驗則是為了與第三組試驗結果對比分析,檢測WACCM6的預報初值在沒有進行同化改善的情況下,對歷史過程的回報能力。

3 2020年北極平流層臭氧3、4月同化模擬與5、6月預報預測

我們首先分析了第一組試驗,也就是WACCM6+DRAT模型在中高層大氣同化了SABER溫度、臭氧和水汽3種變量以后,平流層臭氧的模擬結果。圖1給出了模擬的北極臭氧柱總量(TCO)隨時間的變化特征。圖1a–e分別表示3月1日、3月15日、3月30日、4月15日,4月30日TCO的變化。可見,從3月1日開始,北極圈內已經開始出現臭氧損耗(圖1a);到15日的時候,臭氧損耗繼續加強(圖1b);到了3月30日,北極臭氧已經開始出現大規模損耗現象,TCO降低到了220 DU(Dobson Unit)以下(圖1c);在4月15日時,北極臭氧損耗現象已經開始減弱(圖1d);到4月30日,北極的臭氧大規模損耗事件基本結束(圖1e)。圖2給出了由MERRA2再分析資料計算得到的、整個北極地區平均的TCO隨時間的變化曲線。2020年的臭氧大規模損耗事件主要從2月底開始并到4月中下旬結束(Hu,2020)。3月中下旬,臭氧含量降至最低值,這種創紀錄的低臭氧持續到2020年4月19日。這里需要注意的是,圖2的曲線中TCO的最低值在320 DU左右,這是因為圖2中TCO曲線是60oN~90oN平均得到的。作為參考,在這段時間內,整個北極地區TCO氣候態的平均臭氧總量柱400 DU以上,這是因為北極極渦通常在3月破裂,使富含臭氧的空氣從低緯度地區進入極地地區。自1979年以來,北極上空如此異常的低臭氧只在1997年和2011兩年出現過。2020年3月的異常低臭氧量在觀測記錄中是絕無僅有的。對比圖1和圖2可以看到,WACCM6+DRAT模型的同化試驗模擬結果很好地重現了2020年北極平流層臭氧3、4月大規模損耗情況。

圖1 第一組試驗(同化SABER資料的試驗,見表1)模擬的2020年(a)3月1日、(b)3月15日、(c)3月30日、(d)4月15日和(e)4月30日北極臭氧柱總量(TCO,單位:DU)分布Fig. 1 First set of experiments (Table 1) and simulated total column ozone (TCO, units: DU) distributions for (a) March 1, (b) March 15,(c) March 30, (d) April 15, and (e) April 30, 2020

圖2 MERRA2再分析資料中,過去30年(1990~2020年)60oN~90oN平均的TCO(單位:DU)月變化曲線。黑線為1990~2019年TCO月變化的平均結果,陰影區域代表1990~2019年TCO的月變化范圍,藍線代表2020年TCO月變化曲線Fig. 2 MERRA2 reanalysis of monthly TCO (units: DU) changes averaged over 60°–90° N for the past 30 years (1990–2020). The black line indicates the average of the monthly TCO changes from 1990–2019, while the shaded area represents the range of TCO changes from 1990–2019. The blue line indicates the monthly TCO change for 2020

圖3給出了第二組試驗WACCM6模式在沒有進行資料同化的情況下,模擬的北極臭氧柱總量(TCO)隨時間的變化特征。圖3a–e分別表示3月1日、3月15日、3月30日、4月15日,4月30日TCO的變化。在3月1日和3月15日(圖3a和b),可以看到北極臭氧也出現了損耗的現象。但是,臭氧損耗最嚴重的區域,TCO也在280 DU以上。到3月30日的時候(圖3c),臭氧損耗現象已經完全消失,直到4月底也再沒有出現臭氧損耗現象(圖3d和e)。這說明,沒有同化SABER資料的第二組試驗并沒有重現2020年3、4月北極平流層臭氧大規模損耗的現象。這也說明了,對于目前即使是包含了完整平流層過程的模式模擬的平流層變化與觀測之間還是存在一定的差距。而中高層大氣同化技術可能是解決這一問題的一個有效手段。

圖3 第二組試驗(未進行同化的試驗,見表1)模擬的2020年(a)3月1日、(b)3月15日、(c)3月30日、(d)4月15日和(e)4月30日TCO(單位:DU)分布Fig. 3 Second set of experiments (see Table 1) simulated TCO (units: DU) distributions for (a) March 1, (b) March 15, (c) March 30, (d) April 15,and (e) April 30, 2020

為了進一步的分析歷史重現試驗中(第一、二組試驗)同化試驗和未進行同化的試驗與觀測結果之間的差異,圖4a和圖4b給出了北極地區平均的平流層臭氧含量從2020年3月1日到4月30日之間的變化特征,并與MLS衛星觀測資料進行了對比。對于北極地區,臭氧高含量區域分布在50~1 hPa之間。圖4a、b表明同化試驗模擬的臭氧含量(圖4a)比未進行同化的試驗模擬的臭氧含量(圖4b)要低得多。圖4c顯示了圖4a中的臭氧含量變化與過去20年平均的3月1日到4月30日臭氧含量變化之差。可以看到,相對于過去臭氧的含量變化,2020年3月1日到4月30日北極平流層臭氧的含量要低得多。圖4d顯示了圖4a中的臭氧含量變化與MLS觀測臭氧的差異。總體來說,同化后模擬的平流層臭氧變化與觀測資料結果很接近,在1 hPa以上和100 hPa以下,它們之間的差異非常小。但是,30~1 hPa之間,同化后模擬的臭氧含量相對觀測值也偏小約0.5 ppmv,在100~30 hPa之間,相對觀測值偏大。

圖4 (a)第一組試驗(同化SABER資料的試驗,見表1)模擬的2020年3、4月北極地區平均的臭氧含量變化。(b)第二組試驗(未進行同化的試驗,見表1)模擬的2020年3、4月北極地區平均的臭氧含量變化。(c)為(a)中的臭氧含量變化與過去20年北極地區平均的臭氧變化的氣候態的差值。(d)為(a)中的臭氧含量變化與MLS資料中的北極地區平均的臭氧含量變化的差值。單位:ppm(ppm=10?6)Fig. 4 (a) Simulated changes in Arctic-averaged ozone from the first set of experiments in March and April 2020 (Table 1). (b) Simulated changes in Arctic-averaged ozone from the second set of experiments in March and April 2020 (Table 1). (c) Difference between the ozone changes in (a) and Arctic-averaged ozone change over the past decade. (d) Difference between the ozone changes in (a) and the Arctic-averaged ozone changes based on MLS data. Units: ppm

可見,同化SABER衛星資料有效地提高了WACCM6模擬平流層臭氧變化的能力。在試驗中,一共同化了臭氧、溫度和水汽三種變量。DART的集成卡爾曼濾波器,會將WACCM6的臭氧、溫度和水汽輸出結果調整為與觀測數據信息更一致的狀態。這將可能從兩個方面對WACCM6的臭氧輸出結果進行改善。一是,通過集成卡爾曼濾波器,模擬結果直接向觀測結果的調整。但是,SABER衛星資料軌道覆蓋率不全面,僅僅靠模擬結果直接向觀測結果的調整,對于一些無臭氧觀測資料同化的區域,臭氧模擬改善效果有限。而研究中,溫度和水汽觀測資料的同化會從化學過程方面對臭氧模擬結果進行影響。改善平流層溫度可以改善平流層臭氧的損耗和生成速度模擬。特別對于極地,平流層溫度是影響極地云形成的最重要因子之一。同化平流層溫度資料,使模擬的溫度更接近觀測結果,這也會使極地云的模擬更接近觀測結果。極地云上發生的非均相化學反應釋放的活性鹵族元素,是造成臭氧洞的直接因素。圖1中,第一組同化試驗中能重現2020年春季的臭氧洞現象,除了與直接同化臭氧有重要關系外,同化平流層溫度變化也可能是重要因素。而對于同化平流層水汽資料,水汽分解產生的OH離子可以把臭氧還原成氧氣,另外它也是生成極地云的主要成分。這些都說明,通過同化技術對平流層溫度和水汽模擬的改善,對平流層臭氧模擬的改善作用是肯定的。這里需要說明的是,本研究中,我們僅同化了一種衛星資料。這可能意味著,對于改善平流層變化的模擬,雖然中高層大氣同化技術是一個有效手段,但是僅同化一種衛星資料模擬的平流層變化可能還會與真實大氣之間存在一定誤差。未來還需要同化多種衛星資料,通過多源資料同化的方式來進一步改善平流層模擬效果。

以第一組同化試驗輸出的4月30日分析場為初始場,第三組試驗對2020年5月1日到6月30日的平流層臭氧變化進行了回報試驗檢驗。圖5a給出了回報的平流層臭氧含量變化。圖5b為回報試驗結果與MLS臭氧觀測結果之差。對于0~3天的預報結果,預報值與MLS觀測值之間只有在100~1 hPa存在明顯的差異,在30~1 hPa之間,預測值偏小0.5 ppmv,在100~30 hPa之間,預測值偏大0.5 ppmv;絕對誤差在0.5 ppmv之內。對比圖4d可以看到,差異主要來自于同化試驗輸出的4月30日分析場在100~1 hPa之間與MLS觀測值之間存在誤差。對從第4到第15天的預報結果,預報值與觀測值之間沒有出現誤差隨時間增益的現象,預測值與觀測值之間有很好的一致性。特別是對于第16到第30天的預報結果,其預報值與觀測值之間的差異相比于0~15的預報結果與觀測值之間的差異更小了。對比圖4a和圖5a可以看到,從5月開始,100~1 hPa之間臭氧含量由于季節變化減小。15~30天預報值與觀測值之間差異的減少,應該是由于臭氧總含量隨季節變化減少導致的。也就是說15~30天預報值與觀測值之間的絕對誤差減少了,但是相對誤差可能并沒有顯著減少。從第30天開始的短期氣候尺度預測中的預測值與觀測值之間的誤差突然出現了較大的差異,但是到了45天以后,誤差又減少了(圖5b)。此結果可能與MLS觀測值本身存在的觀測誤差有關。這個結果將來還需要與其他觀測資料對比以進一步驗證。總體來說,短期預測尺度的結果與觀測值之間的誤差很穩定,預測結果與觀測值有很高的一致性。

圖5 (a)第三組試驗(以第一組同化試驗輸出的分析場為初值的試驗,見表1)預報的2020年5月和預測的6月北極地區平均的臭氧含量變化。(b)為(a)圖中的臭氧含量變化與MLS資料中臭氧含量變化的差值。單位:ppmFig. 5 (a) Arctic-averaged ozone changes forecasted and predicted by the third set of experiments (Table 1) for May and June 2020.(b) Difference between ozone changes in (a) and those in the MLS data.Units: ppm

Ivy et al.(2017)和Xie et al.(2016, 2017)發現,北極平流層臭氧損耗能顯著影響北半球氣候變化。利用數值同化技術提升目前模式中模擬、預報和預測平流層臭氧變化,這不但可能提升模式對流層氣候變化的模擬能力,也可能提升對流層延伸期預報和短期氣候預測能力。

4 2020年北極平流層溫度和環流3、4月同化模擬與5、6月預報預測

上一章分析了WACCM6+DRAT同化、預報預測模型對2020年發生的北極臭氧大規模損耗過程的重現能力以及對臭氧的預報預測能力。由于臭氧變化可以通過輻射作用影響平流層的溫度和環流,這一章節,我們將繼續分析WACCM6+DRAT模型對北極平流層溫度和環流2020年3、4月同化模擬和5、6月預報預測的結果。

圖6a、b和c分別顯示了,第一組同化了SABER溫度、臭氧和水汽資料的試驗輸出的、第二組未進行資料同化試驗輸出的以及MERRA2再分析資料的2020年3、4月平流層溫度變化。相比于圖6c,同化試驗和未進行同化的試驗都較好地模擬出了北極2020年3、4月平流層溫度的變化特征(圖6a和b),即平流層低層溫度從3月開始升高,從190 K上升到210 K左右;而平流層中上層溫度變化不大。圖6d顯示了第一組試驗輸出的平流層溫度變化與MERRA2再分析資料之間的差值,圖6e顯示了第二組試驗輸出的平流層溫度變化與MERRA2再分析資料之間的差值。可以看到,在同化了SABER溫度、臭氧和水汽資料以后,在1 hPa以下的溫度變化與MERRA2再分析資料之間的差異小于2 K。但是在1 hPa以上,第一組同化試驗模擬的3月份溫度變化相對于MERRA2再分析資料結果偏大了4 K左右,而4月份卻偏小4 K左右,這可能與WACCM6模式模擬的平流層高層的semi-annual oscillation(SAO)性能有關。但是總體而言,同化后模擬的溫度變化和再分析資料有很高的一致性。未進行同化的試驗模擬的平流層溫度變化與MERRA2再分析資料之間存在較大差異。例如,在3月中上旬,模擬的平流層溫度偏暖6 K左右,而從中下旬開始,平流層下層溫度變化模擬仍然偏暖而上層模擬偏冷 8 K左右。這種差異一直持續到了4月底。

圖6 (a)第一組試驗(同化SABER資料的試驗,見表1)模擬的2020年3、4月北極地區平均的溫度變化。(b)第二組試驗(未進行同化的試驗,見表1)模擬的2020年3、4月北極地區平均的溫度變化。(c)MERRA2再分析資料中的2020年3、4月北極地區平均的溫度變化。(d)為(a)與(c)的差值。(e)為(b)與(c)的差值。單位:KFig. 6 (a) Arctic-averaged temperature change in March and April 2020 simulated in the first set of experiments (Table 1). (b) Arctic-averaged temperature change in March and April 2020 simulated in the second set of experiments (Table 1). (c) MERRA2 reanalysis for Arctic-averaged temperature change in March and April 2020. (d) Difference between (a) and (c). (e) Difference between (b) and (c). Units: K

對應北極平流層溫度的變化,圖7a、b和c分別顯示了,第一組同化了SABER溫度、臭氧和水汽資料的試驗輸出的,第二組未進行資料同化試驗輸出的以及MERRA2再分析資料的2020年3、4月北極平流層環流變化。北極平流層緯向西風在3月初最強,從三月中旬開始減弱,這說明平流層高層極渦從3月中旬開始崩潰,但是平流層低層的極渦還在繼續維持,直到四月中旬以后,緯向西風幾乎接近于零。第一組同化模擬試驗和再分析資料之間表現出了很好的一致性(圖7a和c)。但是,未進行數值同化的第二組試驗模擬得到北極平流層緯向風強度與再分析資料結果存在明顯的差異。圖7b表明,由于未進行資料同化,第二組試驗模擬的整個北極平流層極渦從3月中旬開始就全部崩潰,此后模擬的平流層低層緯向風速比再分析資料里的結果要弱得多。圖7d和圖7e進一步展示了,第一組和第二組試驗模擬的2020年3、4月北極平流層環流變化與再分析資料之間的差異。圖7d所示,第一組同化試驗模擬的1 hPa以下北極平流層環流的變化與再分析資料之間存在高度的一致性;但是,在1 hPa以上,它們之間存在類似日周期的波動性差異。這可能與模式模擬潮汐波性能有一定關系。對于第二組未進行同化的試驗,相對于再分析資料,北極平流層高層的緯向風模擬明顯偏強,但是中低層嚴重偏弱。

圖7 和圖6類似,但是為緯向風變化。單位:m s?1Fig. 7 Same as Fig. 6, but for zonal wind changes. Units: m s?1

同化SABER衛星資料有效地提高了WACCM6模擬平流層溫度和環流變化的能力。同化臭氧、溫度和水汽三種變量,這將可能從兩個方面對WACCM6的溫度輸出結果進行改善。一是,溫度模擬結果直接向溫度觀測結果的調整。臭氧和水汽分別可以吸收短波輻射和長波輻射調整大氣溫度。二則是,臭氧和水汽的模擬通過同化改善后,會通過調整輻射過程改善平流層的溫度模擬。對于平流層的緯向繞極環流,平流層從極地到熱帶溫度梯度是一個主要決定因子。由于平流層熱帶溫度高極地溫度低,平流層溫度從熱帶向極地遞減。根據熱成風原理,平流層會產生繞極環流。上面提到,同化SABER衛星資料提高了WACCM6平流層溫度變化模擬,這也會改善平流層從熱帶到極地溫度梯度的模擬,從而改善平流層緯向風的模擬。另外,對比圖4、圖6和圖7可以發現,由于第二組試驗未同化SABER溫度、臭氧和水汽,其模擬的北極平流層溫度在上層區域偏冷下層區域偏暖,其導致的與中緯度溫度梯度的變化,引起的直接結果是平流層上層緯向風速偏強而中下層偏弱。這對應著北極極渦模擬偏弱。因此,第二組未進行同化的試驗模擬的2020年3、4月北極臭氧損耗比第一組同化試驗和觀測資料顯示的臭氧損耗要弱得多。這再次證明利用同化技術是提升目前模式中平流層模擬能力的一個重要手段。

以第一組同化試驗輸出的4月30日分析場為初始場,第三組試驗對2020年5月1日到6月30日的北極平流層溫度變化進行了回報試驗檢驗。同樣,以第二組未進行同化的試驗輸出的4月30日分析場為初始場,第四組試驗對2020年5月1日到6月30日的平流層溫度變化進行了回報試驗檢驗。為了對比第一組和第二組試驗為第三組和第四組試驗提供的初始場的差別,這里圖8給出了第一組和第二組試驗輸出的以及MERRA2資料中的2020年4月30日北極地區(60°~90°N)平均的溫度和緯向風垂直曲線。對于溫度場來說,在未同化衛星資料的情況下,第二組試驗提供的溫度初始值明顯與MERRA2再分析資料之間存在較大差異,在平流層低層溫度偏小6 k而在高層偏小近8 k(圖8a和b)。而第一組同化了SABER衛星資料溫度的試驗提供的溫度初始值與MERRA2再分析資料之間差異基本在1 hPa以下小于1 K(圖8a和b)。注意1 hPa以上,第一組試驗提供的溫度初始場與MERRA2資料之間也存在2~4 K左右的差異,具體原因將來值得進一步研究。對于緯向風而言,第二組試驗提供的緯向風初始值在平流層低層風速偏小約8 m s?1而在高層偏大約 4 m s?1(圖8c和d)。而第一組試驗提供的緯向風初始值與MERRA2再分析資料之間差異基本在1 hPa以下小于1 m s?1(圖8c和d)。同樣,在1 hPa以上,第一組試驗提供的風速初始場與MERRA2資料之間也存在6 m s?1左右的差異。但是,總體而言,同化試驗提供的初始場相對于未同化的試驗有較大幅度的改善。

圖8 (a)2020年4月30日北極地區(60°~90°N)平均的溫度垂直曲線。黑線基于MERRA2資料,藍線基于第一組試驗(同化SABER資料的試驗)資料,紅線基于第二組試驗(未進行同化的試驗)資料。(b)中藍線為(a)中藍線與黑線之差,代表第一組試驗結果與MERRA2的差值,紅線為(a)中紅線與黑線之差,代表第二組試驗結果與MERRA2的差值。(c)、(d)與(a)、(b)類似,但是為緯向風的變化Fig. 8 (a) Average temperature vertical curve for the Arctic region (60°–90° N) for April 30, 2020. The black line corresponds to MERRA2, the blue line to the first set of experiments, and the red line to the second set of experiments. (b) Difference between the first set of experiments and the MERRA2 (blue line) and difference between the second set of experiments and the MERRA2 (red line). (c) and (d) are similar to (a) and (b), except for zonal wind changes

圖9a和圖9b給出了兩組試驗回報的北極平流層溫度變化,圖9c為MERRA2再分析資料結果。從它們之間的差異可以看到(圖9d和e),以同化試驗輸出的分析場為初始場的平流層溫度預報結果,不管是0~30天的天氣尺度預報還是30~60天的短期氣候尺度預測結果與再分析資料之間在1 hPa以下都非常一致,誤差在2 K以內(圖9d)。在1 hPa以上,0~60天的回報溫度相對再分析資料都偏低,最大偏差大約?6 K左右,特別是4~15天的預報結果相對再分析資料偏差最大。這主要是因為,第一組同化試驗中,模擬的4月底的北極平流層高層溫度相對再分析資料偏低造成的(圖8a和b)。以此作為預測初值,造成了回報試驗中1 hPa以上溫度回報結果相對再分析資料結果偏低(圖9d)。但是有意思的是,此偏差在后續的預測中并沒有被放大而是穩定在了約6 K左右。如前文所述,未進行同化的第二組試驗模擬的4月北極平流層溫度變化與再分析資料之間存在較大的差異,特別是在30 hPa以上,模擬的平流層溫度嚴重偏低(圖8a和b)。因此,以第二組試驗輸出的4月底分析場為初值的第四組試驗,回報的北極平流層溫度在平流層高層和中層都出現了相對再分析資料結果偏低的現象(圖9e)。

圖9 (a)第三組試驗(以第一組同化試驗輸出的分析場為初值的試驗,見表1)預報的2020年5月和預測的6月北極地區平均的溫度變化。(b)第四組試驗(以第二組未進行同化的試驗輸出的分析場為初值的試驗,見表1)預報的5月和預測6月北極地區平均的溫度變化。(c)MERRA2再分析資料中的2020年5、6月北極地區平均的溫度變化。(d)為(a)與(c)的差值。(e)為(b)與(c)的差值。單位:KFig. 9 (a) Forecasted and predicted changes in Arctic-averaged temperature for May and June 2020 by the third set of experiments (Table 1).(b) Forecasted and predicted changes in Arctic-averaged temperature for May and June by the fourth set of experiments (Table 1). (c) Arctic-averaged temperature changes for May and June 2020 based on MERRA2 reanalysis data. (d) Difference between (a) and (c). (e) Difference between (b) and (c).Unis: K

與圖9類似,圖10顯示了北極平流層緯向風變化的回報試驗檢驗。可以看到由于季節變化,從5月份開始,北極平流層高層已經從緯向西風轉變為緯向東風(圖10a、b和c)。對于緯向風的回報,第三組使用同化數據為初值的回報試驗結果與再分析資料之間的差異較小,存在2 m s?1左右的差異(圖10d)。但是在5月初的北極平流層低層,出現了較大幅度的偏差,這可能與第三組試驗對對流層的預報結果存在較大偏差有關。對于第四組沒有使用同化試驗輸出的分析場為初值的試驗,預報預測結果和再分析資料之間存在明顯差異(圖10e)。對于0~15天的天氣尺度預報,北極平流層中上層緯向東風預報偏弱4 m s?1左右,而平流層低層緯向西風預報同樣偏弱4 m s?1左右。這主要與未進行數值同化的試驗(第二組試驗)提供的初值場與再分析資料間在平流層上層和低層存在差異有關(圖7e)。而對于15~30天的預報和30天以上的短期預測結果,預報預測的北極平流層上層緯向東風偏弱(圖10e),差異甚至達到了5 m s?1以上。

圖10 同圖9,但為緯向風變化。單位:m s?1Fig. 10 Same as Fig. 9, but for zonal wind changes. Units: m s?1

5 數值同化技術對全球平流層溫度和環流模擬、預報和預測能力的改善

圖11顯示了,第一組同化了SABER溫度、臭氧和水汽資料的試驗和第二組未進行資料同化試驗輸出的2020年3、4月全球平流層溫度和緯向風與MERRA2再分析資料之間的均方根誤差(RMSE)(圖11a–d)。同 化SABER資 料 后,在1 hPa以下,模擬的全球平流層溫度和緯向風的RMSE分別小于2 K和2 m s?1(圖11a和c)。在1 hPa以上,溫度和緯向風的RMSE分別超過了4 K和6 m s?1。平流層高層與再分析資料之間的差異可能與SABER資料與MERRA2資料之間本身存在差異有一定關系。圖11b和圖11d給出了沒有同化SABER衛星資料的試驗輸出的全球平流層溫度和緯向風的RMSE。可以看到,全球平流層中,模擬的大部分區域的溫度的RMSE超過了4 K,而緯向風的RMSE甚至超過了10 m s?1。對比圖11a和圖11b以及圖11c和圖11d可以看到,同化中高層大氣衛星資料對模擬全球平流層大氣變化都有很顯著的改善效果,而不僅僅局限于上面討論的北極平流層。

圖11 (a)和(c)分別為第一組試驗(同化SABER資料的試驗,見表1)模擬的2020年3、4月平流層溫度和緯向風與MERRA2再分析資料之間的RMSE。(b)和(d)分別為第二組試驗(未進行同化的試驗,見表1)模擬的2020年3、4月平流層溫度和緯向風與MERRA2再分析資料之間的RMSEFig. 11 (a) Stratospheric temperature and (c) wind RMSEs of March and April 2020 between the first set of experiments (Table 1) and the MERRA2 reanalysis data; (b) stratospheric temperature and (d) wind RMSEs of March and April 2020 between the second set of experiments (Table 1) and the MERRA2 reanalysis data

中高空大氣數值同化不但可以改善模式全球平流層大氣變化的模擬能力也可以改善平流層大氣變化的預報與預測性能。圖12顯示了,第三、四組試驗天氣尺度預報和短期氣候尺度預測的2020年5、6月全球平流層溫度與MERRA2再分析資料之間的RMSE。第三組試驗的初值來自于第一組同化試驗輸出的分析場,而第四組試驗的初值來自于第二組未進行同化的試驗輸出的分析場。對0~3天的天氣尺度預報(圖12a和b),第三組試驗預報的平流層溫度的RMSE只在小部分區域超過了4 K,在大部分區域RMSE都小于2 K。而第四組試驗預報的平流層溫度的RMSE在南半球高緯度地區的平流層高層超過了10 K;另外,北極平流層的溫度RMSE也較大,這個特征在前文的圖9也可以看到。第三組試驗天氣尺度預報和短期氣候尺度預測的平流層溫度的RMSE分布相比短期預報并沒有發生變化,雖然RMSE數值有所增大,但是增益較小且穩定(圖12c、e、和g)。而第四組試驗天氣尺度預報和短期氣候尺度預測的平流層溫度的RMSE(圖12d、f、和h)相比第三組試驗的結果要大得多。

圖12 第三組試驗(以第一組同化試驗輸出的分析場為初值的試驗,見表1)輸出的平流層溫度2020年5月(a)0~3天預報、(c)4~15天預報和(e)16~30天預報以及(g)2020年6月短期氣候預測結果與MERRA2再分析資料之間的RMSE。(b)、(d)、(f)和(h)分別為第四組試驗(以第二組未進行同化的試驗輸出的分析場為初值的試驗,見表1)輸出的平流層溫度2020年5月(b)0~3天預報、(d)4~15天預報和(f)16~30天預報以及(h)6月短期氣候預測結果與MERRA2再分析資料之間的RMSEFig. 12 Stratospheric temperature RMSEs between the (a) 0–3 days forecast, (c) 4–15 days forecast, (e) 16–30 days forecast, and (g) short-term climate prediction from the third set of experiments (Table 1) and the MERRA2 reanalysis data; stratospheric temperature RMSEs between the (b) 0–3 days forecast, (d) 4–15 days forecast, (f) 16–30 days forecast, and (h) short-term climate prediction from the fourth set of experiments (Table 1) and the MERRA2 reanalysis data

圖13顯示了,第三、四組試驗天氣尺度預報和短期氣候尺度預測的2020年5、6月全球平流層緯向風與MERRA2再分析資料之間的RMSE。第三組試驗0~3天的天氣尺度預報的平流層風速的RMSE較小,15天以后的預報和短期氣候預測的緯向風的RMSE有所增大;特別在熱帶和南半球平流層高層,緯向風的預報預測結果與MERRA2再分析資料之間存在較大差異;北半球平流層環流預測結果相對較好(圖13a、c、e和g)。而第四組試驗全球平流層緯向風預報和預測結果的RMSE都很大(圖13b、d、f和h)。總體來說,使用同化試驗輸出的分析場作為初值,能較好地改善模式平流層溫度和環流的預報預測性能。

圖13 同圖12,但為緯向風RMSEFig. 13 Same as Fig. 12, but for zonal wind RMSE

6 討論與結論

作為臨近空間的重要組成部分,平流層是未來空間飛行器的主要活動區域。平流層的環境大氣狀態對飛行器的準確入軌和安全都具有顯著影響,也是飛行器的設計參數、飛行試驗的主要依據。但是,目前平流層的觀測資料相對較少,而數值模式在表達平流層環境大氣狀態時還存在一定缺陷。數據同化技術可將新的觀測數據引入模型,這有利于模擬減少甚至濾掉模型噪聲,使得模擬預測結果更加接近平流層環境真實狀態。同化技術與模型的結合,可以獲得許多沒有觀測儀器支持下的大氣狀態數據。因此,平流層的數據同化技術成為了連接觀測數據和模型模擬預測的主要橋梁,在平流層的狀態表征與研究中會起到重要作用。平流層數據同化、預報預測系統的研究,將對平流層數值預報質量的提高起到關鍵作用,從而加強對平流層低動態飛行器飛行活動的氣象安全保障。本研究使用WACCM6模式+DART工具,開發了中高層大氣溫度、臭氧和水汽資料的同化接口,搭建了包含完整平流層過程的數值同化、天氣預報和短期氣候預測模型。利用該模型,針對2020年3、4月北極平流層臭氧大規模損耗事件進行了重現模擬試驗,并以同化試驗輸出的分析場作為初值,對5~6月的平流層大氣進行了0~30天天氣尺度預報以及31~60天短期氣候尺度預測。結果表明:本模型能真實反映2020年3、4月北極平流層出現的大規模臭氧損耗事件隨時間的演變特征,和MLS衛星觀測結果很接近;利用同化試驗輸出的4月末分析場作為初值,預報的5月北極平流層臭氧變化與MLS衛星觀測值的差值小于0.5 ppmv,預測的6月北極平流層臭氧變化只有在30~10 hPa之間區域,與觀測之間的差異超過了1 ppmv。該模型同化模擬的3~4月、預報預測的5~6月北極平流層溫度和緯向風變化與MERRA2再分析資料結果具有很好的一致性,僅在北極平流層頂部,預報預測的溫度和緯向風與再分析資料之間的RMSE分別為約3 K和約4 m s?1。與未進行數值同化的試驗相比,中高層大氣同化技術對平流層中低層模擬效果改善最為顯著,其預報預測結果比未進行同化的試驗的預報預測結果的誤差減少了50%以上。

目前該研究還存在以下一些問題沒有討論:

(1)文章只描述了該模型對平流層環流等變量同化、預報和預測效果的好與壞,并沒有解釋導致這些效果的原因,這也是下一步工作需要詳細研究的問題;

(2)由于臭氧和水汽有很強的輻射作用,它們的變化是影響平流層溫度變化的主要因子;平流層臭氧和水汽的變化又和溫度有關。本研究同時同化了溫度、臭氧、水汽3種變量,從而不能清楚判定哪一個變量的同化對平流層臭氧、溫度和環流模擬的改善效果最顯著;

(3)本研究只同化了一種衛星資料,多種衛星資料同時同化對模式的改善作用也有待進一步研究;

(4)我們目前在模擬時間上只針對性地選取2020年3~6月進行了研究,本模型對不同月份不同季節的改善效果還不清楚。這些問題都需要進一步的討論和研究。

致謝感謝NCAR提供的WACCM6模式和DART同化工具;感謝NASA提供的MERRA2、MLS和SABER資料;感謝北京師范大學超算中心提供的計算資源。

猜你喜歡
大氣
大氣的呵護
軍事文摘(2023年10期)2023-06-09 09:15:06
首次發現系外行星大氣中存在CO2
科學(2022年5期)2022-12-29 09:48:56
宏偉大氣,氣勢與細膩兼備 Vivid Audio Giya G3 S2
太赫茲大氣臨邊探測儀遙感中高層大氣風仿真
有“心氣”才大氣
如何“看清”大氣中的二氧化碳
學生天地(2020年18期)2020-08-25 09:29:24
大氣穩健的美式之風Polk Audio Signature系列
稚拙率真 圓融大氣
中國篆刻(2017年3期)2017-05-17 06:20:46
大氣古樸揮灑自如
大氣、水之后,土十條來了
新農業(2016年18期)2016-08-16 03:28:27
主站蜘蛛池模板: 97精品久久久大香线焦| 美女一级免费毛片| 亚洲中文精品人人永久免费| 狠狠干综合| 亚洲永久免费网站| 99久久人妻精品免费二区| 午夜福利亚洲精品| 欧美在线黄| 欧美日本二区| 亚洲综合中文字幕国产精品欧美 | 97久久免费视频| 久久精品亚洲中文字幕乱码| 六月婷婷激情综合| 毛片久久久| 欧美成人一级| 国产最新无码专区在线| 日本精品影院| 日韩大片免费观看视频播放| 久久亚洲中文字幕精品一区| 国产在线精品人成导航| 看av免费毛片手机播放| 日本欧美视频在线观看| 澳门av无码| 中文字幕 91| 国产69囗曝护士吞精在线视频| 美女无遮挡免费视频网站| 日韩国产一区二区三区无码| 国产成人亚洲精品蜜芽影院| 伊人激情综合| 亚洲精品777| 欧美a在线| 精品一区二区三区视频免费观看| 无码专区第一页| 成人日韩视频| 国产午夜人做人免费视频中文| 午夜无码一区二区三区| 日韩一区二区三免费高清| 欧美a在线视频| 99re视频在线| 99精品免费在线| 亚洲开心婷婷中文字幕| 日本伊人色综合网| 国产91小视频| 9啪在线视频| 野花国产精品入口| 久久频这里精品99香蕉久网址| 欧美a级完整在线观看| 亚洲成年人网| 伊人久久久大香线蕉综合直播| 国产毛片不卡| 亚洲无码高清免费视频亚洲 | 国内毛片视频| 日韩AV无码一区| 成人91在线| 久久国产高清视频| 91久久性奴调教国产免费| 日韩精品免费在线视频| 秋霞一区二区三区| 亚洲Va中文字幕久久一区| 97se亚洲综合在线| 波多野结衣无码视频在线观看| 欧美区一区二区三| 国产成人综合久久| 亚洲精品欧美重口| 狠狠色婷婷丁香综合久久韩国| 亚洲系列无码专区偷窥无码| 97久久人人超碰国产精品| 日本午夜视频在线观看| 99在线国产| 成人毛片免费观看| 欧美视频在线第一页| 韩国福利一区| 国产丝袜啪啪| 性喷潮久久久久久久久| 亚洲精品无码日韩国产不卡| 婷婷综合缴情亚洲五月伊| 天天色综网| 性欧美在线| 国产制服丝袜无码视频| 人妻一区二区三区无码精品一区| 亚洲免费播放| 久久精品国产999大香线焦|