亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        兩次典型沙塵的WRF-Chem數(shù)值模擬:不同起沙方案模擬效果的評(píng)估

        2022-11-08 06:12:00譚成好賈世國唐明金
        地球化學(xué) 2022年5期

        尹 馨, 譚成好, 賈世國, 唐明金

        兩次典型沙塵的WRF-Chem數(shù)值模擬:不同起沙方案模擬效果的評(píng)估

        尹 馨1, 3, 譚成好1, 3, 賈世國2*, 唐明金1

        (1. 中國科學(xué)院 廣州地球化學(xué)研究所, 有機(jī)地球化學(xué)國家重點(diǎn)實(shí)驗(yàn)室/廣東省環(huán)境資源利用與保護(hù)重點(diǎn)實(shí)驗(yàn)室/粵港澳大灣區(qū)環(huán)境污染過程與控制聯(lián)合實(shí)驗(yàn)室, 廣東 廣州 510640; 2. 中山大學(xué) 大氣科學(xué)學(xué)院, 廣東 廣州 510275; 3. 中國科學(xué)院大學(xué), 北京 100049)

        為探尋更適用于中國沙塵過程數(shù)值模擬的起沙方案, 利用WRF-Chem模式分別耦合GOCART、AFWA和Shao04 3種不同起沙參數(shù)化方案對(duì)2017年5月3–5日和2018年3月26–29日發(fā)生在我國西北地區(qū)的2次典型沙塵過程進(jìn)行數(shù)值模擬, 并綜合多種地面和衛(wèi)星觀測(cè)資料對(duì)模擬結(jié)果進(jìn)行評(píng)估。結(jié)果表明: 通過3種起沙方案計(jì)算得到的起沙量差距較大, 其中GOCART方案在0~6 μm粒徑的起沙量最大、Shao04方案最小。從近地面PM10濃度來看, 3種方案均能反映出2次事件中典型城市PM10濃度的變化特征, Shao04方案模擬的PM10濃度和觀測(cè)相關(guān)性最高; 從氣溶膠光學(xué)厚度(AOD)來看, 通過與AERONET、MODIS和CALIPSO的觀測(cè)資料綜合比較, 發(fā)現(xiàn)模擬結(jié)果對(duì)AOD存在系統(tǒng)性低估, GOCART方案與觀測(cè)結(jié)果最近似, 其次為Shao04方案, 最后為AFWA方案。綜合考慮認(rèn)為, GOCART起沙方案相對(duì)更適用于我國西北地區(qū)的沙塵模擬。

        沙塵; WRF-Chem; 起沙參數(shù)化方案

        0 引 言

        沙塵顆粒物作為全球主要的大氣氣溶膠成分之一, 年排放量可達(dá)約5000 Tg(Tegen and Fung, 1994)。沙塵天氣不僅對(duì)人類健康、空氣質(zhì)量和氣候變化產(chǎn)生直接作用, 還會(huì)通過散射、吸收太陽短波輻射與地面長波輻射和作為云凝結(jié)核參與成云致雨等方式直接或間接地影響地氣系統(tǒng)的輻射收支平衡(劉曉東等, 2004; Shao et al., 2011; Sun et al., 2012; Choobari et al., 2014; Huang et al., 2014; 王訓(xùn)明等., 2015)。此外, 當(dāng)攜帶有大量生物可利用關(guān)鍵微量元素的沙塵氣溶膠被遠(yuǎn)距離輸送到海洋中時(shí)(Rea et al., 1998, Shao et al., 2011), 沙塵顆粒物還會(huì)影響生物地球化學(xué)循環(huán)過程(Jickells et al., 2005; Mahowald et al., 2006; 高會(huì)旺等, 2009; Maher et al., 2010)。

        鑒于沙塵循環(huán)在氣候系統(tǒng)中發(fā)揮的重要作用, 國內(nèi)外學(xué)者曾采用地基觀測(cè)、衛(wèi)星遙感、統(tǒng)計(jì)分析、數(shù)值模擬等手段對(duì)沙塵進(jìn)行了系列研究(Huebert et al., 2003; Mikami et al., 2006; Huang et al., 2008; Che et al., 2015; Xia et al., 2016)。近年來, 數(shù)值模式愈發(fā)成為研究和預(yù)報(bào)沙塵天氣的有力工具, 被廣泛應(yīng)用于研究沙塵的排放、傳輸、沉降及其對(duì)環(huán)境和氣候的影響。當(dāng)前數(shù)值模式可以較好地還原氣象條件和沙塵基本特征, 然而不同起沙參數(shù)化方案模擬的起沙量具有顯著差異, 主要原因在于不同方案對(duì)起沙過程模擬的差異大。不同的起沙參數(shù)化方案考慮的輸入?yún)?shù)、氣象和土壤參數(shù)基準(zhǔn)不一致、計(jì)算簡化方式不同, 因而不同起沙方案模擬的起沙量和沙塵濃度的時(shí)空分布等存在較大的不確定性(Huneeus et al., 2010; Sugimoto et al., 2013; 吳成來和林朝暉, 2014; Kim and Choi, 2015; Su and Fung, 2015; Flaounas et al., 2017, 康麗泰和陳思宇, 2017, 劉筱冉等, 2018, Yuan et al., 2019)。此外, 不同起沙方案在不同地區(qū)的適應(yīng)性存在一定的差異。西北干旱和半干旱地區(qū)是我國沙塵天氣的高發(fā)區(qū), 也是東亞沙塵的主要來源, 沙塵顆粒物經(jīng)遠(yuǎn)距離傳輸可對(duì)位于華北、東北和長江中下游的城市空氣質(zhì)量產(chǎn)生影響, 然而當(dāng)前對(duì)中國西北干旱、半干旱區(qū)域的沙塵數(shù)值模擬研究中針對(duì)不同起沙方案的評(píng)估對(duì)比有限, 因此探究起沙方案對(duì)于中國西北沙塵事件模擬的影響具有重要意義。

        本研究選取2017年5月3–5日和2018年3月26–29日期間發(fā)生在我國西北地區(qū)的2次典型沙塵天氣過程, 使用WRF-Chem模式分別耦合在該模式中最常用的3種起沙參數(shù)化方案GOCART、AFWA和Shao04進(jìn)行數(shù)值模擬, 并綜合多種地面和衛(wèi)星觀測(cè)資料進(jìn)行對(duì)比評(píng)估。目標(biāo)是探討不同起沙方案模擬的異同, 為優(yōu)化模式性能、深入理解和探究沙塵循環(huán)的大氣物理化學(xué)過程提供科學(xué)依據(jù)。

        1 數(shù)據(jù)與方法

        1.1 WRF-Chem簡介

        WRF(weather research and forecasting model)模式系統(tǒng)是由美國氣象界聯(lián)合開發(fā)的新一代中尺度預(yù)報(bào)模式和同化系統(tǒng), 常被用于全球的區(qū)域尺度、中小尺度數(shù)值模擬、物理參數(shù)化方案研究、實(shí)時(shí)業(yè)務(wù)天氣預(yù)報(bào)、數(shù)據(jù)同化研究、區(qū)域氣候模擬以及耦合空氣質(zhì)量模式。WRF-Chem模式是將WRF與化學(xué)模塊(Chemistry, 簡稱Chem)實(shí)現(xiàn)完全在線耦合的區(qū)域大氣動(dòng)力–化學(xué)模式(Grell et al., 2005), 其在線同步計(jì)算物理和化學(xué)過程, 考慮了生物體排放、人類活動(dòng)排放和氣相化學(xué)與氣溶膠化學(xué)、輸送(包括平流、擴(kuò)散和對(duì)流過程)、干濕沉降等過程,能真實(shí)、同步地反映大氣物理化學(xué)過程。因此, WRF-Chem模式不僅能模擬污染過程中的風(fēng)速、氣溫和氣壓等氣象要素的特征, 還能模擬大氣污染物的排放、傳輸和分布等, 是用于沙塵模擬研究的有力工具。

        1.2 起沙參數(shù)化方案

        1.2.1 GOCART方案

        GOCART(goddard chemistry aerosol radiation and transport)方案(Ginoux et al., 2001)是基于10 m風(fēng)速和土壤濕度的經(jīng)驗(yàn)公式, 該方案分5個(gè)粒徑段計(jì)算垂直起沙通量, 分別為0~1、1~1.8、1.8~3、3~6和6~10 μm, 沙塵粒子的垂直起沙通量計(jì)算公式如下:

        式中:為垂直起沙通量;為常量0.8 μg·s?2·m5;為風(fēng)蝕指數(shù), 代表潛在沙塵源分布受地表特征如冰雪覆蓋率和植被率等的影響;s為各粒徑段沙塵的質(zhì)量比例, 除0~1 μm粒徑為0.1, 其他各粒徑均為0.25;10為離地10 m的水平風(fēng)速;t為臨界摩擦速度, 受沙塵粒徑大小和土壤濕度影響。

        1.2.2 AFWA方案

        AFWA(air force weather agency)方案(Jones et al., 2010, 2012)是GOCART方案的改進(jìn)方案,簡稱為AFWA方案, 繼承并發(fā)展了Marticorena and Bergametti (MB)公式中(Marticorena and Bergametti, 1995)關(guān)于跳沙轟擊的起沙機(jī)制計(jì)算方法:

        式中:(D)為躍移水平沙通量;為可變動(dòng)的經(jīng)驗(yàn)比例常量, 本研究中值為1;ρ為空氣密度;為重力加速度;*為摩擦速度;*t為臨界摩擦速度;是所有粒徑的總轟炸量;bulk是垂直沙塵通量;為風(fēng)蝕指數(shù);為考慮了黏粒比例的沙塵轟擊率。首先根據(jù)式(2)和式(3)計(jì)算總起沙通量, 再基于脆性材料斷裂理論將沙塵分類到5個(gè)粒徑段(0~1、1~1.8、1.8~3、3~6和6~10 μm, 質(zhì)量比分別為0.1074、0.1012、0.2078、0.4817和0.1019)(Kok, 2011), 從而可計(jì)算出各粒徑的垂直起沙通量。和GOCART方案相比, AFWA方案使用摩擦速度*和臨界摩擦速度*t取代了10 m風(fēng)速, 當(dāng)土壤濕度高于干燥臨界點(diǎn)時(shí)增大, 反之不作調(diào)整。也就是說, AFWA方案是在GOCART方案的基礎(chǔ)上對(duì)黏粒比例、風(fēng)蝕指數(shù)和臨界摩擦速度等主要參數(shù)進(jìn)行了改進(jìn)。

        1.2.3 Shao04方案

        Shao04方案(Shao, 2004)屬于UoC(the university of cologne)方案中的1種, 該方案兼顧了跳沙轟擊和轟擊粒子分裂兩種起沙機(jī)制, 方案計(jì)算結(jié)果與觀測(cè)的一致性較好。Shao04方案計(jì)算過程如下:

        式中:(d,d)為由粒徑為d的沙粒跳躍碰撞產(chǎn)生的粒徑為d的垂直起沙通量;c為比例系數(shù);η,為由土壤排放到大氣中的沙粒的質(zhì)量分?jǐn)?shù); σ是轟炸系數(shù);為與沙塵粒徑分布有關(guān)的權(quán)重因子, 表示集合沙粒被釋放的難易程度;Q為躍移水平沙塵通量;是所有粒徑的總垂直起沙通量;是第種粒徑分布的占比;σ和均用于第種粒徑的對(duì)數(shù)分布;ρ是土壤容重;p()和p()分別為粒徑的最小和全分布, 二者被用于計(jì)算土壤的塑性壓力p()。

        因此, Shao04方案計(jì)算各粒徑沙塵通量的準(zhǔn)確性主要依賴于原始土壤數(shù)據(jù)和。和AFWA方案相似, Shao04方案同樣使用了沙暴參數(shù)化來計(jì)算垂直起沙通量, 但計(jì)算過程和粒徑分布(0~2、2~3.6、3.6~6、6~12和12~20 μm 5個(gè)粒徑)和AFWA有所差異。在Shao04方案中, 首先使用式(11)將最大粒徑可達(dá)2000 μm的粒徑分為由小到大的100種粒徑分布, 然后根據(jù)式(6)~(9)計(jì)算各粒徑的垂直起沙通量。隨后, 規(guī)定分類直徑為2、3.6、6、12和20 μm, 由此獲得0~2、2~3.6、3.6~6、6~12和12~20 μm的垂直起沙通量。與AFWA和GOCART方案不同, Shao04方案只使用風(fēng)蝕指數(shù)來規(guī)定潛在的沙塵源。

        1.3 數(shù)值模擬實(shí)驗(yàn)設(shè)計(jì)

        2017年5月3–7日, 我國內(nèi)蒙古部分地區(qū)發(fā)生沙塵暴, 局地出現(xiàn)強(qiáng)沙塵暴, 西北地區(qū)、華北地區(qū)和華東多地受到沙塵天氣影響。2018年3月26–29日, 我國西北地區(qū)、華北地區(qū)、華中北部和東北多部先后出現(xiàn)揚(yáng)沙或浮塵, 其中內(nèi)蒙古錫林郭勒盟局地出現(xiàn)沙塵暴。這2次沙塵事件的特征是持續(xù)時(shí)間短、強(qiáng)度大、影響范圍廣, 屬于近年來較為嚴(yán)重的沙塵污染事件。

        本研究運(yùn)用WRF-Chem V3.8.1模擬上述兩次沙塵案例, 模式設(shè)置采用單層模擬, 范圍覆蓋整個(gè)中國地區(qū), 水平分辨率為27 km, 區(qū)域格點(diǎn)數(shù)為250×170,垂直分辨率為38層, 時(shí)間積分步長為150 s。模擬時(shí)間分別為2017年5月1–9日和2018年3月22–31日。模式在化學(xué)模塊方案的選擇上采用耦合了RACM (regional atmospheric chemistry mechanism)化學(xué)反應(yīng)機(jī)制的GOCART氣溶膠方案(Hewson et al., 2012), 同時(shí)不考慮光化學(xué)過程、生物氣溶膠、海鹽氣溶膠、濕清除和氣溶膠對(duì)輻射的反饋等機(jī)制, 主要選用的其他物理化學(xué)方案見表1。

        1.4 數(shù)據(jù)資料

        研究選取的氣象驅(qū)動(dòng)資料為歐洲氣象中心提供的每6 h一次、水平分辨率為1°×1°的ERA-Interim再分析數(shù)據(jù)。人為源采用清華大學(xué)的多尺度排放清單數(shù)據(jù)庫MEIC(multi-resolution emission inventory for China)清單(http://www.meicmodel.org)。為了對(duì)模擬結(jié)果進(jìn)行驗(yàn)證, 本研究采用了酒泉、包頭、北京和青島4個(gè)城市監(jiān)測(cè)點(diǎn)的PM10濃度監(jiān)測(cè)數(shù)據(jù)驗(yàn)證PM10地面濃度模擬、采用AERONET經(jīng)過云處理和人工檢查的2.0級(jí)數(shù)據(jù)驗(yàn)證站點(diǎn)氣溶膠光學(xué)厚度(aerosol optical depth, AOD)模擬, 采用搭載在Aqua衛(wèi)星上MODIS提供的結(jié)合了暗像元和深藍(lán)算法反演的L3級(jí)550 nm每日平均數(shù)據(jù)產(chǎn)品驗(yàn)證AOD的水平分布特征模擬, 采用CALIPSO的Level-2氣溶膠廓線產(chǎn)品驗(yàn)證氣溶膠的垂直分布特征模擬。

        表1 模式主要選用的參數(shù)化方案

        2 結(jié)果與分析

        2.1 起沙量對(duì)比

        前48 h為模式spin-up, 本研究通過歷史天氣記錄驗(yàn)證, 僅分析記錄有沙塵天氣的連續(xù)數(shù)日起沙量, 時(shí)間段分別為2017年5月3–5日和2018年3月27–29日。

        由于3種方案均包含對(duì)0~6 μm粒徑的模擬, 討論該粒徑范圍的垂直起沙量可從同一粒徑水平分析不同方案的起沙特征, 表2可見兩次案例在0~6 μm和全粒徑的模擬起沙量, 起沙量的計(jì)算包括了模擬區(qū)域內(nèi)全部格點(diǎn)。對(duì)于2017年5月的沙塵案例, 3種方案模擬的起沙量均為從3日到5日逐日減小, 但方案間每日的結(jié)果差異明顯, GOCART方案起沙量逐日降低的趨勢(shì)較平緩、AFWA方案和Shao04方案的起沙量從4日到5日驟降。分析起沙量結(jié)果, GOCART方案的0~6 μm粒徑總起沙量最大(91.64 Tg),其次為AFWA方案(86.64 Tg), Shao04方案在該粒徑的總起沙量顯著低于前兩者(43.53 Tg)。其次, 將GOCART方案和AFWA方案的0~10 μm全粒徑起沙量與Shao04方案0~12 μm起沙量進(jìn)行對(duì)比, 結(jié)果仍為GOCART方案最大(121.95 Tg)、AFWA方案次之(86.64 Tg)、Shao04方案最小(82.19 Tg)。GOCART方案中0~6 μm粒徑的起沙量為全粒徑的75%、AFWA方案為90%, Shao04方案0~6 μm粒徑的起沙量為0~12μm粒徑的53%、僅為0~20 μm全粒徑的33%, 這說明GOCART方案和AFWA方案模擬結(jié)果中0~6 μm粒徑的更小顆粒沙塵所占比例更大。

        表2 2次案例的模擬起沙量(Tg)

        2018年3月沙塵事件, 模擬起沙量的逐日變化趨勢(shì)相同, 均為逐日增大, 但各方案模擬的總起沙量相差較大。GOCART方案在0~6 μm的起沙量模擬值最高(90.18 Tg), AFWA方案約為GOCART方案的30%~50%, Shao04方案最小, 其中在27日的模擬結(jié)果甚至僅為GOCART方案的1%。GOCART方案在0~6 μm的總起沙量是0~10 μm全粒徑總起沙量的75%; AFWA方案在0~6 μm的總起沙量為38.89 Tg, 為0~10 μm粒徑的90%; Shao04方案在0~6 μm的總起沙量為10.29 Tg, 是0~12 μm粒徑的54%、0~20 μm粒徑的35%。對(duì)于2018年3月案例的模擬, 不論是同一粒徑范圍還是單個(gè)方案考慮的全粒徑范圍, Shao04方案的起沙量均為最小。

        綜合2017和2018年兩次沙塵事件的模擬結(jié)果, GOCART、AFWA和Shao04方案間的差異顯著, 因此, 有必要進(jìn)行進(jìn)一步分析, 通過對(duì)比觀測(cè)和模擬結(jié)果來探討不同方案對(duì)我國沙塵模擬的適用性。

        2.2 地面觀測(cè)與模擬對(duì)比

        2.2.1 城市PM10濃度

        為了更全面地評(píng)估不同起沙方案的模擬效果, 本研究選取了鄰近塔克拉瑪干沙漠的酒泉、鄰近戈壁沙漠的包頭、位于下風(fēng)向的北京和沿海傳輸路徑的青島4個(gè)城市的PM10觀測(cè)值來與同時(shí)段的模擬結(jié)果進(jìn)行對(duì)比, 每日觀測(cè)結(jié)果和模擬結(jié)果均進(jìn)行了6 h, 以便更明顯地觀察變化趨勢(shì)。圖1、2分別為典型城市站點(diǎn)2017年5月3–9日和3月25–31日期間模擬與觀測(cè)PM10的對(duì)比。

        2017年5月3–9日期間, 酒泉3日中午地面PM10濃度陡升至1150 μg·m?3, Shao04方案反映了這一高值特征但估值明顯偏低; GOCART方案雖然可以反映高值但是出現(xiàn)時(shí)間向前偏移; AFWA高值不明顯。此后, 觀測(cè)濃度在5日達(dá)到第2個(gè)峰值(400 μg·m?3), 模擬結(jié)果中該峰值均早于觀測(cè)結(jié)果6 h,GOCART和AFWA方案的濃度值與觀測(cè)相近, Shao04方案顯著偏低。對(duì)于酒泉地面PM10濃度的模擬, Shao04方案相較于觀測(cè)結(jié)果的相關(guān)系數(shù)最高(0.89); 3個(gè)方案對(duì)于除峰值以外的其他時(shí)段模擬表現(xiàn)比較相近且基本符合觀測(cè)。

        在包頭, 前12 h的模擬結(jié)果數(shù)值均顯著偏高。觀測(cè)峰值共出現(xiàn)兩次, 分別為1350 μg·m?3(4日06時(shí))和1200 μg·m?3(5日06時(shí)), 兩次峰值的間隔時(shí)間比較短, 3個(gè)方案均模擬還原出了兩個(gè)峰值, 但模擬峰值分別提前了6 h和12 h。6–9日期間, 包頭的觀測(cè)濃度在每天正午達(dá)到小峰值, 但模擬結(jié)果沒有還原出這一特征, 模擬結(jié)果在非峰值時(shí)段偏低。整體來看, 不同起沙方案的PM10模擬結(jié)果在酒泉和包頭等近沙塵源地的差異性并不顯著, 各方案均能較好地模擬出觀測(cè)到的沙塵濃度變化特征, 這說明了模擬是可信的。

        圖1 典型城市站點(diǎn)2017年5月3–9日模擬與觀測(cè)PM10對(duì)比(Obs代表觀測(cè)值, 下同)

        圖2 典型城市站點(diǎn)2018年3月25–31日模擬與觀測(cè)PM10對(duì)比

        北京的觀測(cè)濃度在4–5日間出現(xiàn)較大起伏, PM10濃度變化超過1000 μg·m?3, 其余時(shí)段濃度在300 μg·m?3以下, 和包頭相似。北京的觀測(cè)濃度在較短時(shí)間間隔中也出現(xiàn)了2次峰值(4日12時(shí)和5日06時(shí), 均為1000 μg·m?3左右)。GOCART和AFWA方案在北京的差異性小, 二者的統(tǒng)計(jì)分析數(shù)據(jù)非常接近(與觀測(cè)數(shù)據(jù)的相關(guān)系數(shù)分別為0.48、0.46), Shao04方案對(duì)兩次峰值的模擬濃度均高于GOCART和AFWA方案, 對(duì)第1個(gè)峰值的模擬出現(xiàn)一定程度的估值偏高但更接近觀測(cè)值, 模擬結(jié)果的第2次峰值, 3種方案均估值偏低, 其中Shao04方案與觀測(cè)值差值最小, 但Shao04方案在模擬后期出現(xiàn)1次虛高值。

        青島的觀測(cè)濃度了出現(xiàn)兩次時(shí)間間隔較長的峰值, 分別在6日0時(shí)和7日12時(shí), PM10濃度同約為400 μg·m?3。GOCART和AFWA方案在青島的差異性不顯著, 3個(gè)起沙方案均捕捉到了首次峰值, 其中GOCART方案和AFWA方案的模擬濃度值近似但提前12 h出現(xiàn), Shao04方案僅提前6 h但濃度上約為觀測(cè)值的2倍, 模擬結(jié)果均沒有表現(xiàn)出第2次峰值。

        2018年3月25–31日期間, 酒泉的觀測(cè)濃度有3次顯著高值出現(xiàn), 分別為900(26日0時(shí))、1000(29日18時(shí))和2300 μg·m?3(31日0時(shí)), 其余時(shí)間段濃度值在400 μg·m?3以下。模擬結(jié)果均沒有捕捉到首次峰值, 28日12時(shí)前模擬濃度值小于200 μg·m?3。后兩次峰值方案間的差異性較大: GOCART方案分別提前18 h、6 h出現(xiàn)峰值, 且2次峰值濃度都在800 μg·m?3以下; AFWA方案均提前6 h, 2次濃度均低于500 μg·m?3; Shao04方案分別提前6 h、12 h, 濃度分別為500 μg·m?3、2000 μ·m?3, 該方案對(duì)峰值的模擬更接近觀測(cè)值。

        包頭的PM10觀測(cè)濃度在28日12時(shí)前持續(xù)增大, 從起始的約150 μg·m?3至該時(shí)刻出現(xiàn)首次峰值約900 μg·m?3, 第2次峰值于29日18時(shí)出現(xiàn)(520 μg·m?3),最后1次峰值于31日0時(shí)出現(xiàn)(375 μg·m?3)。GOCART方案捕捉到了首次和末次峰值, 但模擬出的首次峰值提前了36 h而末次峰值的模擬值(700 μg·m?3)比觀測(cè)值高約1.5倍; AFWA方案的模擬濃度值始終不超過300 μg·m?3, 對(duì)峰值的捕捉效果不佳; Shao04方案模擬出了3次完整峰值, 但峰值時(shí)間均提前(分別提前12、18和18 h), 且模擬值明顯偏低。綜合分析酒泉和包頭兩地的PM10濃度, GOCART方案下的模擬結(jié)果整體偏高, 而AFWA方案和Shao04方案偏低, 這和起沙量模擬結(jié)果中GOCART方案顯著高于AFWA和Shao04方案的現(xiàn)象一致。

        北京的觀測(cè)濃度呈現(xiàn)單峰趨勢(shì), 峰值于28日12時(shí)出現(xiàn)(1600 μg·m?3), 其余時(shí)段濃度值較平穩(wěn)(不超過400 μg·m?3)。3個(gè)方案的模擬結(jié)果在28日18時(shí)均出現(xiàn)了升高趨勢(shì), 但濃度和觀測(cè)值相差甚遠(yuǎn), 其中Shao04方案在29日18時(shí)出現(xiàn)了和觀測(cè)結(jié)果不符的峰值(900 μg·m?3)??傮w來看, 模擬低估了北京的地面PM10濃度。

        青島的PM10觀測(cè)濃度始終小于200 μg·m?3, 說明該次事件中青島受沙塵影響不顯著, 雖然濃度變化范圍較大但是特征值不明顯。起沙方案間的差異性在青島無法得以體現(xiàn), 模擬結(jié)果的濃度變化范圍大于觀測(cè)值, 但整體相差不大。

        總的來說, 3種方案均能反映出2次事件中典型城市PM10濃度的變化特征, 但各方案對(duì)兩次沙塵的模擬表現(xiàn)出了不同特征。對(duì)于2017年沙塵, 方案間的PM10模擬結(jié)果在近沙塵源地的兩個(gè)西北城市差異性并不顯著, 在下風(fēng)向城市Shao04方案存在估值偏高; 對(duì)于2018年沙塵, GOCART方案對(duì)2個(gè)西北城市的PM10濃度的估值偏高, 3個(gè)方案對(duì)下風(fēng)向城市的模擬差異不顯著。另外, Shao04方案的2次模擬結(jié)果在北京均出現(xiàn)1個(gè)沙塵后期虛高峰值。

        2.2.2 AERONET站點(diǎn)觀測(cè)

        AOD反映了整個(gè)大氣層或某特定大氣層中顆粒物對(duì)太陽輻射的削弱程度。2017年5月沙塵期間, Beijing-CAMS站點(diǎn)(116.3°E, 39.9°N, 106.0 m)含5個(gè)AOD有效數(shù)據(jù)日, 4日和8日數(shù)據(jù)缺失(圖3a)。3日的AOD最大(1.6), 而模擬AOD值均在0.2以下, 模擬值在4日出現(xiàn)峰值(Shao04方案0.9、GOCART方案0.6、AFWA方案0.5), 后降至較低的平穩(wěn)值(不超過0.2), 7~9日的增大趨勢(shì)和觀測(cè)相符; XiangHe站點(diǎn)(117.0°E, 39.8°N, 36.0 m), 含5個(gè)有效數(shù)據(jù)日, 5日和8日數(shù)據(jù)缺失(圖3b)。觀測(cè)AOD值在4日最大, 為2.2、7日最小, 為0.1, 基本呈單峰趨勢(shì)。模擬準(zhǔn)確還原了觀測(cè)的單峰趨勢(shì), 但相較于觀測(cè)結(jié)果數(shù)值顯著偏低, 最大值僅為1.0。綜合2個(gè)AERONET站點(diǎn)結(jié)果, 模擬基本反映了AOD的變化趨勢(shì)。

        2018年3月沙塵期間, Beijing-CAMS站點(diǎn)含6個(gè)有效數(shù)據(jù)日, 28日缺乏數(shù)據(jù), 觀測(cè)AOD值范圍在0.2~1.4之間, 模擬結(jié)果在30日前存在普遍偏低(模擬AOD值不超過0.3), 同時(shí)模擬沒有反映出此次沙塵期間站點(diǎn)AOD趨勢(shì)的變化, 27日觀測(cè)AOD值一度高達(dá)1.4但模擬未能還原這一高值; 31日, 觀測(cè)AOD值由0.6增大至0.9, 模擬均反映了這一增大趨勢(shì)(從0.2增大至超過0.8)??傮w來說, 對(duì)于2018年沙塵, 各方案間AOD模擬差異不明顯, 對(duì)該次沙塵的AOD反映效果較差(圖4)。

        圖3 2017年5月3–9日AERONET觀測(cè)AOD(550 nm波段)與模擬對(duì)比

        2.3 衛(wèi)星觀測(cè)與模擬對(duì)比

        2.3.1 MODIS衛(wèi)星觀測(cè)

        MODIS衛(wèi)星觀測(cè)數(shù)據(jù)可反映沙塵的水平空間分布特征, 圖5為2017年5月4日MODIS觀測(cè)和WRF-Chem模擬的550 nm 波段AOD水平分布對(duì)比。內(nèi)蒙古、華北和東北地區(qū)在4日受沙塵影響嚴(yán)重, AOD值超過1.0(圖5a)。對(duì)比觀測(cè)和模擬的AOD水平分布結(jié)果, 發(fā)現(xiàn)3種起沙方案基本均反映了AOD水平分布特征, 但方案間有所不同。3種方案均在華北地區(qū)出現(xiàn)高值區(qū), 但GOCART模擬的高值中心往西南偏移, AFWA方案高值區(qū)域面積較小, Shao04方案模擬的高值區(qū)最接近觀測(cè), 但高值范圍較觀測(cè)偏小, 未模擬出甘肅、內(nèi)蒙古北部的次高值區(qū), 模擬AOD水平分布所反映的特征和各方案的起沙量模擬數(shù)值有較好的對(duì)應(yīng)關(guān)系, 方案間起沙量的模擬差異是造成AOD的水平分布差異的原因之一??傮w上看, 模擬結(jié)果較觀測(cè)值偏低, 特別是在內(nèi)蒙古, 模擬值約為0.4, 而觀測(cè)值在1.0以上。

        圖4 2018年3月25–31日AERONET觀測(cè)AOD(550 nm波段)與模擬對(duì)比

        圖6為2018年3月29日MODIS觀測(cè)和WRF-Chem模擬的550 nm 波段AOD水平分布對(duì)比。29日, 我國有2個(gè)明顯的AOD高值區(qū), 分別位于新疆的塔克拉瑪干沙漠和山東半島, 可知沙塵在起沙源地的影響仍然在持續(xù)。而前期起沙已擴(kuò)散至海上, 對(duì)朝鮮半島和日本產(chǎn)生了影響, 這一特征在模擬結(jié)果中均有體現(xiàn)。對(duì)比觀測(cè)和模擬結(jié)果, GOCART方案和AFWA方案均體現(xiàn)了塔克拉瑪干沙漠的高值特征, 其中AFWA方案模擬AOD值較低, 而Shao04方案沒有模擬出這一區(qū)域的高值。3個(gè)方案在戈壁沙漠均有較高值出現(xiàn), 但缺少觀測(cè)結(jié)果無從對(duì)比。觀測(cè)結(jié)果中存在跨河南–山東半島–朝鮮半島的帶狀A(yù)OD高值區(qū)域(超過1.0), 模擬結(jié)果中反映了這一特征, 但模擬值僅約0.4, 和2017年結(jié)果相同, 模擬AOD值較觀測(cè)值偏低。

        圖5 2017年5月4日MODIS觀測(cè)AOD(550 nm波段)與模擬對(duì)比(紅色直線為該事件中對(duì)比的CALIPSO衛(wèi)星軌跡)

        圖6 2018年3月29日MODIS觀測(cè)AOD(550 nm波段)與模擬對(duì)比(紅色直線為該事件中對(duì)比的CALIPSO衛(wèi)星軌跡)

        2.3.2 CALIPSO衛(wèi)星觀測(cè)

        CALIPSO衛(wèi)星數(shù)據(jù)可提供氣溶膠的垂直結(jié)構(gòu)信息, 圖7a反映了2017年5月3日06:00–06:09期間CALIPSO衛(wèi)星觀測(cè)到的532 nm消光系數(shù)。受模式輸出步長限制, 選取與該時(shí)段最接近的06:00時(shí)刻模式輸出結(jié)果, 不難發(fā)現(xiàn)模擬基本還原了氣溶膠的垂直結(jié)構(gòu)特征: CALIPSO觀測(cè)和WRF-Chem模擬的最大擴(kuò)散高度均在8.0 km左右; 觀測(cè)結(jié)果在38°~46°N垂直截面的消光系數(shù)出現(xiàn)高值區(qū), 擴(kuò)散高度最大達(dá)到約4.5 km。模擬結(jié)果捕捉到了這一高值區(qū), 最大擴(kuò)散高度略低于觀測(cè)值(4.0 km), 其中Shao04方案的模擬結(jié)果和觀測(cè)到的垂直結(jié)構(gòu)相似度最高; 模擬方案在32~36°N有明顯的垂直混合特征, 其中GOCART的消光系數(shù)達(dá)到0.3以上, 同區(qū)域缺少觀測(cè)值。

        圖8為2018年3月28日觀測(cè)(19:33–19:45)與模擬(20:00)的氣溶膠垂直結(jié)構(gòu)對(duì)比。CALIPSO觀測(cè)中出現(xiàn)40~43°N的高值區(qū), 模擬反映了這一垂直擴(kuò)散特征, 但缺少地形和觀測(cè)值, 限制了可對(duì)比范圍。GOCART方案下的模擬結(jié)果和觀測(cè)最為近似, 二者的相關(guān)系數(shù)最大值均大于0.4、擴(kuò)散高度均超過3.0 km;AFWA方案模擬的消光系數(shù)最小、擴(kuò)散高度在2.5 km以下; Shao04方案的擴(kuò)散趨勢(shì)和AFWA方案相似, 但消光系數(shù)還原度高于AFWA方案。整體來說, GOCART方案表現(xiàn)最好。

        3 結(jié) 論

        (1) 對(duì)于2次沙塵案例, 3種起沙方案下的起沙量結(jié)果差異較大, 其中GOCART方案在0~6 μm粒徑的起沙量最大、Shao04方案最小。相較于GOCART方案和AFWA方案, Shao04方案對(duì)2018年沙塵的起沙量模擬明顯偏小。

        (2) 3種方案均能反映出2次事件中典型城市PM10的變化特征。對(duì)于近沙塵源地的PM10模擬, 各方案對(duì)2017年沙塵事件的模擬差異不顯著, 而GOCART方案對(duì)2018年案例的模擬值明顯偏高。整體來看, Shao04方案模擬的PM10濃度和觀測(cè)相關(guān)性最高。

        (3) 通過與AERONET、MODIS和CALIPSO 3種觀測(cè)資料的對(duì)比分析, 發(fā)現(xiàn)模擬結(jié)果對(duì)AOD存在系統(tǒng)性估值偏低, 同時(shí), 模擬AOD水平分布所反映的特征與各方案的起沙量模擬數(shù)值有較好的對(duì)應(yīng)關(guān)系, 模擬結(jié)果中方案間起沙量的差異是造成AOD的水平分布差異的原因之一。Shao04方案對(duì)于塔克拉瑪干沙漠的起沙量模擬的數(shù)值偏低且在塔克拉瑪干沙漠區(qū)域的沙塵水平分布特征模擬效果不顯著。GOCART方案對(duì)2次案例特征的整體模擬效果最好, 其次為Shao04方案, 最后為AFWA方案。

        圖7 2017年5月3日CALIPSO觀測(cè)與模擬的氣溶膠垂直結(jié)構(gòu)(橫坐標(biāo)表示觀測(cè)點(diǎn)的經(jīng)緯度)

        圖8 2018年3月28日CALIPSO觀測(cè)與模擬的氣溶膠垂直結(jié)構(gòu)

        (4) 綜合2次事件模擬結(jié)果與觀測(cè)數(shù)據(jù)的對(duì)比分析, 發(fā)現(xiàn)GOCART方案和Shao04方案均優(yōu)于AFWA方案且具備各自的優(yōu)勢(shì), 但Shao04方案在我國最主要的沙塵源地——塔克拉瑪干沙漠的起沙模擬不確定性較大, 因此認(rèn)為GOCART起沙方案相對(duì)更適用于我國西北地區(qū)的沙塵模擬。

        致謝:特別感謝南京大學(xué)黃昕副教授和另外一名匿名審稿專家對(duì)本文提出的建設(shè)性的修改意見及建議。

        高會(huì)旺, 祁建華, 石金輝, 石廣玉, 馮士筰. 2009. 亞洲沙塵的遠(yuǎn)距離輸送及對(duì)海洋生態(tài)系統(tǒng)的影響. 地球科學(xué)進(jìn)展, 24(1): 1–10.

        康麗泰, 陳思宇. 2017. 中國北方一次沙塵天氣過程的數(shù)值模擬. 中國沙漠, 37(2): 321–331.

        劉曉東, 田良, 張小曳. 2004. 塔克拉瑪干沙塵活動(dòng)對(duì)下游大氣PM10濃度的影響. 中國環(huán)境科學(xué), 24(5): 17–21.

        劉筱冉, 王金艷, 邱繼勇, 李全喜, 魏林波. 2018. 起沙方案對(duì)西北地區(qū)沙塵過程模擬的影響. 環(huán)境保護(hù)科學(xué), 44(4): 69–76.

        王訓(xùn)明, 周娜, 郎麗麗, 花婷, 焦琳琳, 馬文勇. 2015. 風(fēng)沙活動(dòng)對(duì)陸地生態(tài)系統(tǒng)影響研究進(jìn)展. 地球科學(xué)進(jìn)展, 30(6): 627–635.

        吳成來, 林朝暉. 2014. WRF/Chem模式中兩種起沙參數(shù)化方案對(duì)東亞地區(qū)一次強(qiáng)沙塵暴過程模擬的影響. 氣候與環(huán)境研究, 19(4): 419–436.

        Che H, Zhang X Y, Xia X, Goloub P, Holben B, Zhao H, Wang Y, Zhang X C, Wang H, Blarel L, Damiri B, Zhang R, Deng X, Ma Y, Wang T, Geng F, Qi B, Zhu J, Yu J, Shi G. 2015. Ground-based aerosol climatology of China: Aerosol optical depths from the China aerosol remote sensing network (CARSNET) 2002–2013., 15(13): 7619–7652.

        Chen S H, Sun W Y. 2002. A one-dimensional time dependent cloud model., 80(1): 99–118.

        Chen Y, Yang K, Zhou D, Qin J, Guo X. 2010. Improving the Noah land surface model in arid regions with an appropriate parameterization of the thermal roughness length., 11: 995–1006.

        Choobari O A, Zawar-Reza P, Sturman A. 2014. The global distribution of mineral dust and its impacts on the climatesystem: A review., 138: 152–165.

        Flaounas E, Kotroni V, Lagouvardos K, Klose M, Flamant C, Giannaros T M. 2017. Sensitivity of the WRF-Chem (V3.6.1) model to different dust emission parametrisation: Assessment in the broader Mediterranean region., 10(8): 2925–2945.

        Ginoux P, Chin M, Tegen I, Prospero J, Holben B, Dubovik O, Lin S J. 2001. Sources and distributions of dust aerosols simulated with the GOCART model., 106: 20255–20274.

        Grell G A, Peckham S E, Schmitz R, McKeen S A, Frost G, Skamarock W C, Eder B. 2005. Fully coupled “online” chemistry within the WRF model., 39(37): 6957–6975.

        Grell G A, Dévényi D. 2002. A generalized approach to parameterizing convection combining ensemble and data assimilation techniques., 29(6): 587–590.

        Hewson M, McGowan H, Phinn S. 2012. Comparing remotely sensed and modelled aerosol properties for a region of low aerosol optical depth. 2012 IEEE International Geoscience and Remote Sensing Symposium. Munich, IGARSS: 2512–2515.

        Hong S Y, Noh Y, Dudhia J. 2006. A new vertical diffusion package with an explicit treatment of entrainment processes., 134(9): 2318–2341.

        Huang J, Minnis P, Chen B, Huang Z, Liu Z, Zhao Q, Yi H, Ayers J. 2008. Long-range transport and vertical structure of Asian dust from CALIPSO and surface measurements during PACDEX., 113, D23212.

        Huang J, Wang T, Wang W, Li Z, Yan H. 2014. Climate effects of dust aerosols over East Asian arid and semiarid regions., 119(19): 11398–11416.

        Huebert B J, Bates T, Russell P B, Shi G Y, Kim Y L, Kawamura K, Carmichael G, Nakajima T. 2003. An overview of ACE-Asia: Strategies for quantifying the relationships between Asian aerosols and their climatic impacts., 108, D23.

        Huneeus N, Michael S, Balkanski Y, J G, S K, Prospero J, S B, Boucher O, Chin M, Dentener F, Diehl T, Easter R, D F, Ghan S, Ginoux P, A G, Horowitz L, Koch D, Krol M, Zender C. 20110. Global dust model intercomparison in AEROCOM., 11: 7781–7816.

        Jickells T D, An Z S, Andersen K K, Baker A R, Bergametti G, Brooks N, Cao J J, Boyd P W, Duce R A, Hunter K A, Kawahata H, Kubilay N, Laroche J, Liss P S, Mahowald N, Prospero J M, Ridgwell A J, Tegen I, Torres R. 2005. Global iron connections between desert dust, ocean biogeochemistry, and climate., 308(5718): 67– 71.

        Jones S L, Adams-Selin R, Hunt E D, Creighton G A, Cetola J D. 2012. Update on modifications to WRF-CHEM GOCART for fine-scale dust forecasting at AFWA // American Geophysical Union Fall Meeting Abstracts. Washington D C: AGU Publication, A33D-0188.

        Jones S L, Creighton G A, Kuchera E L, George K D, Elliott A J. 2010. Adapting WRF-CHEM GOCART for fine-scale dust forecasting // American Geophysical Union Fall Meeting Abstracts. Washington D C: AGU Publication, U14A-06.

        Kim H L, Choi M H. 2015. Impact of soil moisture on dust outbreaks in East Asia: Using satellite and assimilation data., 42(8): 2789–2796.

        Kim Y, Fu J S, Miller T L. 2010. Improving ozone modeling in complex terrain at a fine grid resolution: Part I — examination of analysis nudging and all PBL schemes associated with LSMs in meteorological model., 44(4): 523–532.

        Kok J F. 2011. A scaling theory for the size distribution of emitted dust aerosols suggests climate models underestimate the size of the global dust cycle., 108(3): 1016–1021.

        Lin Y L, Farley R D, Orville H D. 1983. Bulk parameterization of the snow field in a cloud model., 22(6): 1065–1092.

        Maher B A, Prospero J M, Mackie D, Gaiero D, Hesse P P, Balkanski Y. 2010. Global connections between aeolian dust, climate and ocean biogeochemistry at the present day and at the last glacial maximum., 99(1–2): 61–97.

        Mahowald N M, Muhs D R, Levis S, Rasch P J, Yoshioka M, Zender C S, Luo C. 2006. Change in atmospheric mineral aerosols in response to climate: Last glacial period, preindustrial, modern, and doubled carbon dioxide climates., 111, D10202.

        Marticorena B, Bergametti G. 1995. Modeling the atmospheric dust cycle: 1. Design of a soil-derived dust emission scheme., 100(D8): 16415–16430.

        Mielikainen J, Huang B, Huang A, Goldberg M D. 2012. GPU acceleration of the updated goddard shortwave radiation scheme in the weather research and forecasting (WRF) model., 5(2): 555–562.

        Mikami M, Shi G Y, Uno I, Yabuki S, Iwasaka Y, Yasui M, Aoki T, Tanaka T Y, Kurosaki Y, Masuda K, Uchiyama A, Matsuki A, Sakai T, Takemi T, Nakawo M, Seino N, Ishizuka M, Satake S, Fujita K, Hara Y, Suzuki J. 2006. Aeolian dust experiment on climate impact: An overviewof Japan-China joint project ADEC., 52(1–4): 142–172.

        Rea D K, Snoeckx H, Joseph L H. 1998. Late Cenozoic eolian deposition in the North Pacific: Asian drying, Tibetan uplift, and cooling of the Northern Hemisphere., 13(3): 215–224.

        Shao Y. 2004. Simplification of dust emission scheme and comparison with data., 109, D10202.

        Shao Y P, Wyrwoll K H, Chappell A, Huang J P, Lin Z H, McTainsh G H, Mikami M, Tanaka T Y, Wang X L, Yoon S C. 2011. Dust cycle: An emerging core theme in Earth system science., 2(4): 181–204.

        Su L, Fung J C H. 2015. Sensitivities of WRF-Chem to dust emission schemes and land surface properties in simulating dust cycles during springtime over East Asia., 120(21): 11215–11230.

        Sugimoto N, Hara Y, Shimizu A, Nishizawa T, Matsui I, Nishikawa M. 2013. Analysis of dust events in 2008 and 2009 using the lidar network, surface observations and the CFORS model., 49: 27–39.

        Sun Y, Clemens S C, Morrill C, Lin X P, Wang X L, An Z S. 2012. Influence of Atlantic meridional overturning circulation on the East Asian winter monsoon., 5(1): 46–49.

        Tegen I, Fung I. 1994. Modeling of mineral dust in the atmosphere: Sources, transport, and optical thickness., 99(D11): 22897–22914.

        Xia X, Che H, Zhu J, Chen H, Cong Z, Deng X, Fan X, Fu Y, Goloub P, Jiang H, Liu Q, Mai B, Wang P, Wu Y, Zhang J, Zhang R, Zhang X. 2016. Ground-based remote sensing of aerosol climatology in China: Aerosol opticalproperties, direct radiative effect and its parameterization., 124: 243–251.

        Yuan T G, Chen S Y, Huang J P, Zhang X R, Luo Y, Ma X J, Zhang G L. 2019. Sensitivity of simulating a dust storm over Central Asia to different dust schemes using the WRF-Chem model., 207: 16–29.

        Simulating two typical dust storms with the WRF-Chem model:Sensitivity of different dust emission schemes

        YIN Xin1, 3, TAN Chenghao1, 3, JIA Shiguo2*, TANG Mingjin1

        (1. State Key Laboratory of Organic Geochemistry / Guangdong Key Laboratory of Environmental Protection and Resources Utilization / Guangdong-Hong Kong-Macao Joint Laboratory for Environment Pollution and Control,Guangzhou Institute of Geochemistry, Chinese Academy of Sciences, Guangzhou 510640, Guangdong, China; 2. School of Atmospheric Sciences, Sun Yat-Sen University, Guangzhou 510275, Guangdong, China; 3. University of Chinese Academy of Sciences, Beijing 100049, China)

        Dust emission events can cause large fluctuations in particulate matter concentrations and need to be considered to predict atmospheric events accurately. This study evaluated the impact of dust emission schemes on the performance of the WRF-Chem model by comparing the data collected from two severe dust events recorded during 3–5 May 2017 and 26–29 March 2018 over northwest China against three different simulation schemes; specifically, the WRF-Chem model coupled with the GOCART, AFWA and, Shao04 schemes. Results show that there is a significant difference between the simulated results of different dust emission schemes for both events: The GOCART scheme always presented the highest dust emission in the 0–6 μm size range while the Shao04 scheme had the lowest. While the trend of surface PM10concentration in typical cities was accurately predicted by all the dust emission schemes, the Shao04 scheme provided the best performance. Generally, the 3 dust emission schemes underestimated the aerosol optical depth (AOD) when compared with AERONET, MODIS and CALIPSO data. The GOCART scheme performed best, followed by Shao04 and AFWA schemes in AOD simulations. Overall, the GOCART scheme is the best choice for dust simulation in northwest China.

        dust storm; WRF-Chem; dust emission parameterization scheme

        X823; X831

        A

        0379-1726(2022)05-0528-12

        10.19700/j.0379-1726.2022.05.003

        2020-11-03;

        2020-12-11

        國家重點(diǎn)研發(fā)計(jì)劃(2018YFC0213901)和國家自然科學(xué)基金(42022050)聯(lián)合資助。

        尹馨(1996–), 女, 碩士研究生, 環(huán)境工程專業(yè)。E-mail: yinxin18@mails.ucas.ac.cn

        賈世國(1985–), 男, 副教授, 主要從事大氣化學(xué)研究。E-mail: jiashg3@mail.sysu.edu.cn

        波多野结衣免费一区视频| 国产精品自拍网站在线| 日本一道本加勒比东京热| 熟女中文字幕一区二区三区| 各种少妇正面着bbw撒尿视频| 色婷婷五月综合亚洲小说| 天堂最新在线官网av| 亚洲综合国产精品一区二区| 久久精品国产免费观看三人同眠| 亚洲avav天堂av在线网爱情| 爽妇网国产精品| 国产一区二区三区白浆在线观看| 亚洲视频在线观看一区二区三区| 好男人社区影院www| AV无码最在线播放| 中文字幕麻豆一区二区| 午夜国产精品视频在线观看| 免费高清av一区二区三区| 福利一区视频| 亚洲视频精品一区二区三区 | 97在线观看| 中文字幕一区二区三区在线不卡 | 欧美午夜一区二区福利视频| 精品久久久久久国产潘金莲| 亚洲一区二区三区日韩在线观看 | 久久精品国产亚洲av成人文字| 久久婷婷人人澡人人喊人人爽| 久久免费国产精品| 人妻熟女中文字幕在线视频| 亚洲综合av一区二区三区蜜桃| 亚洲日本在线电影| 精品一二区| 91九色熟女潮喷露脸合集| 亚洲乱亚洲乱妇| 水蜜桃久久| 美腿丝袜网址亚洲av| 精品免费国产一区二区三区四区| 欧美黑人粗暴多交高潮水最多| 国产亚洲AV片a区二区| 国产激情视频在线观看大全| 亚洲精品无码国产|