• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    邊界層方案對南京地區(qū)PM2.5濃度模擬的影響

    2021-08-09 02:12:50王安庭李煜斌杜秋燕王曉東高志球南京信息工程大學(xué)大氣物理學(xué)院江蘇南京0044中國科學(xué)技術(shù)大學(xué)地球與空間科學(xué)學(xué)院安徽合肥3006
    中國環(huán)境科學(xué) 2021年7期
    關(guān)鍵詞:邊界層擴散系數(shù)湍流

    王安庭,李煜斌*,趙 純,杜秋燕,王曉東,高志球 (.南京信息工程大學(xué)大氣物理學(xué)院,江蘇 南京 0044;.中國科學(xué)技術(shù)大學(xué)地球與空間科學(xué)學(xué)院,安徽 合肥 3006)

    PM2.5是大氣污染的重要污染物之一[1-3].準(zhǔn)確模擬 PM2.5濃度是開展氣象預(yù)報業(yè)務(wù)、保障生產(chǎn)生活的重要基礎(chǔ)[4-6].大氣污染物主要集中于大氣邊界層[7-10],大氣邊界層中的湍流輸運是污染物擴散的主要機制[11-14].湍流的發(fā)生和發(fā)展受制于大氣邊界層的動力和熱力學(xué)結(jié)構(gòu)[15],如小風(fēng)穩(wěn)定層結(jié)大氣將抑制湍流,進(jìn)而可能導(dǎo)致強污染事件[16-17],而在大風(fēng)或?qū)α餍蕴鞖獾目刂葡?湍流發(fā)展旺盛,通??諝馕廴镜某潭纫蚕鄬^低[18].因此,數(shù)值模式中污染物濃度的模擬結(jié)果對邊界層物理狀態(tài)和物理過程的刻畫十分敏感[19-21].

    國內(nèi)學(xué)者對不同邊界層方案在不同地區(qū)的污染物濃度模擬表現(xiàn)開展了一系列研究.徐敬等[22]和王繼康等[23]在華北地區(qū)對比了 YSU[24]、MYJ[25]和ACM2[26]對臭氧和 PM2.5模擬的影響,均發(fā)現(xiàn)局地閉合的 MYJ方案模擬所得白天邊界層污染物濃度明顯高于具有非局地閉合項的YSU方案和ACM2方案.這一現(xiàn)象在董春卿等[27]于山西地區(qū)對比YSU和MYJ方案的研究中并不明顯,這可能與局地下墊面條件、氣象背景及污染物平流輸送強度有關(guān).蔡子穎等[28]對比了 YSU、MYJ、BL[29]和 MYNN3[30]方案對天津地區(qū)PM2.5的預(yù)報效果,發(fā)現(xiàn)晴空和大風(fēng)天氣BL方案表現(xiàn)較優(yōu),陰天和小風(fēng)天氣則YSU和MYJ方案較為理想.Zhao等[31]對比了 YSU、MYJ和ACM2對京津冀地區(qū)2015年夏季模擬的臭氧濃度,發(fā)現(xiàn)不同邊界層方案可導(dǎo)致超過 20%的濃度模擬差別.

    國外,Cuchiara等[32]采用 YSU、MYJ、ACM2和QNSE[33]方案對德克薩斯州一次污染天氣狀況進(jìn)行模擬,發(fā)現(xiàn)YSU方案對污染物的垂直混合更接近觀測值.Banks等[34]對比YSU、MYJ、ACM2和BL方案對西班牙伊比利亞半島春季的空氣質(zhì)量模擬,發(fā)現(xiàn)ACM2模擬臭氧濃度更為準(zhǔn)確.Parra[35]分析了YSU、MYJ、ACM2、MYNN2[36]、MYNN3和GFS方案[37]對厄瓜多爾安第斯山脈地區(qū)污染物濃度模擬的差別,發(fā)現(xiàn) YSU和 MYJ具有更準(zhǔn)確的模擬能力.De lange等[38]對比了 YSU、MYJ、ACM2和MYNN2在非洲南部高地草原地區(qū)的氣象要素模擬能力,發(fā)現(xiàn)MYJ和MYNN2的模擬結(jié)果更為接近觀測值.

    綜上,不同邊界層方案在不同地區(qū)的表現(xiàn)不一,針對特定地區(qū)深入研究不同邊界層方案的模擬特性十分必要[39-41].依據(jù)局地閉合和非局地閉合可將當(dāng)前邊界層方案分為 3類,即:非局地閉合方案(如YSU)、局地閉合方案(如 MYJ、BL、QNSE、MYNN2、MYNN3)、局地與非局地混合型閉合方案(如ACM2),本文選取此上述中的4種代表性方案YSU、MYJ、MYNN2和 ACM2進(jìn)行了對比.同時,根據(jù)湍流閉合的階數(shù),將其分為:一階閉合(如 YSU 和ACM2)和高階閉合(如1.5階的BL、QNSE、MYJ、MYNN2和MYNN3)兩類.本文選取YSU和ACM2作為1階閉合代表、而MYJ和MYNN2作為高階閉合方案的代表進(jìn)行分析.本文使用這 4種邊界層方案對南京地區(qū)冬季一次典型污染過程進(jìn)行模擬,將模擬結(jié)果與觀測和再分析資料進(jìn)行對比分析,解析不同方案在該地區(qū)的模擬特性及相應(yīng)機制,旨在為該地區(qū)空氣質(zhì)量模擬和預(yù)報提供參考.

    1 資料與方法

    選取2016年12月28日~2017年1月5日之間南京冬季一次較為典型污染過程進(jìn)行研究.期間,南京地區(qū)受一弱高壓中心影響,風(fēng)速較小,盛行西南風(fēng),出現(xiàn)了多次中度污染天氣,且于2016年12月31日中午出現(xiàn)了短暫的重度污染天氣,至2017年1月5日空氣質(zhì)量轉(zhuǎn)為良好.PM2.5濃度則由2016年12月28日的 25μg/m3逐漸增加到 12月 31日的接近200μg/m3,而后回落至 75μg/m3附近.

    1.1 觀測資料

    本文使用的觀測資料包括地面站點資料、探空廓線資料和 ERA5再分析資料.地面站點資料包括時間分辨率為 1h的氣象要素及 PM2.5濃度觀測,來源于“國家環(huán)境空氣質(zhì)量監(jiān)測網(wǎng)”于南京市布置的9個站點,即中華門(118.78°E,32.02°N)、玄武湖(118.81°E,32.07°N)、仙林大學(xué)城(118.92°E,32.09°N)、山西路(118.77°E,32.07°N)、瑞金路(118.81°E,32°04″N)、浦口(118.57°E,32.06°N)、邁皋橋(118.82°E,32.11°N)、草場門(118.75°E,32.06°N)和奧體中心(118.73°E,32.02°N).在分析過程中,本文取9個站點的觀測平均值與臨近格點的模擬數(shù)據(jù)平均值進(jìn)行比較.探空數(shù)據(jù)來源于中國氣象局南京站點(118°48″E,32°00″N)的 L 波段雷達(dá),其每日北京時間8:00和20:00獲取氣象要素的垂直廓線信息,秒探空數(shù)據(jù)垂直分辨率約為 4m(包含溫度和濕度廓線),分探空數(shù)據(jù)約為 300m(包含風(fēng)速風(fēng)向廓線).邊界層高度觀測數(shù)據(jù)來自于一臺型號為LB110-ESS-D200的激光雷達(dá),它發(fā)射355nm波長的激光脈沖,重復(fù)頻率為20hz,在無人看管下可24h工作[42].ERA5再分析資料來源于歐洲中期天氣預(yù)報中心網(wǎng)站(https://cds.climate.copernicus.eu/#!/home),時間分辨率為逐小時,水平分辨率為0.25°×0.25°,垂直分為37層.

    1.2 WRF-Chem參數(shù)設(shè)置

    WRF-Chem模式在美國國家大氣研究中心(NCAR)研發(fā)的天氣預(yù)報模式WRF基礎(chǔ)上加入了化學(xué)模塊,并在時間和空間上將化學(xué)傳輸模塊與氣象模塊完全耦合.其不僅可用于天氣預(yù)報或氣候預(yù)測,還可模擬顆粒物和臭氧等污染物的濃度、擴散和相互作用等[43-44].WRF-Chem模式是研究區(qū)域大氣污染狀況的重要工具[45].本文所采用的WRF-Chem版本是由中國科技大學(xué)更新的WRF-Chem3.5.1版本.與WRF-Chem官網(wǎng)公開發(fā)布的版本相比,中國科技大學(xué)的WRF-Chem版本包含一些附加功能,如氣溶膠-雪相互作用、氣溶膠物種的輻射強迫診斷以及陸面耦合生物源的VOC排放等.

    模式主要設(shè)置如下:模式模擬的區(qū)域如圖 1所示,空間分辨率為 15km,具有 119×112個單元網(wǎng)格.其氣象邊界由1°×1°的NCEP 全球再分析資料FNL提供;針對化學(xué)邊界,本文依據(jù) Du等[46]及其前序論文[47-51],設(shè)計了一個 lat-lon投影、水平分辨率為1°×1°的全球模擬,為本文的模擬區(qū)域提供化學(xué)邊界.垂直方向分層為41層,其中2km以下有18層,最底層約為64m,模式層頂為100hpa.排放源的層數(shù)為20層,氣象化學(xué)機制采用 CBMZ[52-53]方案,光解過程由FAST-J[54]方案在線計算,氣溶膠參數(shù)化選用MOSAIC模型.本文模式設(shè)置是中國科學(xué)技術(shù)大學(xué)WRF-Chem版本模式設(shè)置的延續(xù)[46-50].

    圖1 模式的模擬區(qū)域Fig.1 The simulation area

    模式中微物理過程參數(shù)化方案采用 Morrison two-moment[37]方案,長波輻射和短波輻射方案均采用 RRTMG[55]方案,積云對流參數(shù)化方案采用Kain-Fritsch[56]方案.陸面過程參數(shù)方案采用 Noah方案[57-60].

    本文采用 4種邊界層方案,即 YSU、MYJ、MYNN2和ACM2方案,并增加2個敏感性實驗方案 ACM2R2和 ACM2R5.其中 ACM2R2在 ACM2方案的基礎(chǔ)上,將模式底部6層(約500m)的湍流擴散系數(shù)最小閾值設(shè)定為 2m2/s.而 ACM2R5在ACM2方案的基礎(chǔ)上,將模式底部10層(約1300m)的湍流擴散系數(shù)最小閾值設(shè)定為5m2/s.另外,YSU、MYNN2、ACM2、ACM2R2和ACM2R5近地面層均對應(yīng)MM5表面層方案,而MYJ對應(yīng)Eta表面層方案.

    YSU方案是非局地閉合K理論方案.在湍流局地擴散的基礎(chǔ)上增加了對流性大尺度湍渦導(dǎo)致的非局地混合作用,同時考慮邊界層頂?shù)膴A卷過程.

    MYJ方案是高階局地閉合方案,通過湍流動能來閉合湍流方程,求解邊界層內(nèi)的湍流通量.MYJ考慮的物理過程復(fù)雜,不能滿足假設(shè)條件時模擬結(jié)果有一定偏差,如在某些區(qū)域容易低估大氣邊界層高度[61].

    MYNN2方案具有預(yù)報次網(wǎng)格尺度湍流動能的優(yōu)勢.MYNN2方案將氣壓相關(guān)項進(jìn)行了參數(shù)化處理,同時考慮浮力的影響,使得其邊界層高度模擬具有顯著改善[61].

    ACM2方案特點是通過調(diào)節(jié)湍流擴散項和非局地項之間的比例系數(shù),實現(xiàn)由穩(wěn)定條件下的局地擴散算法到不穩(wěn)定條件下的局地擴散項加非局地擴散項算法的平緩過渡.

    1.3 模式氣象初始場資料及排放源清單

    模式氣象初始場資料采用分辨率為 1°×1°的NCEP 全球再分析資料FNL,其時間分辨率為6h.模擬時間段為北京時間2016年12月28日00:00~2017年1月5日00:00.

    模式外部采用準(zhǔn)全球模擬的人為排放是從第 2版半球污染傳輸(HTAPv2)[62]中獲得的,水平分辨率為 0.1°×0.1°,時間為 2010 年.模式內(nèi)部采用的人為排放源清單來自于清華大學(xué)開發(fā)的中國多尺度排放清單模型 MEIC[63-65]2015 版本,分辨率為 0.1°×0.1°.MEIC清單是一套基于云計算平臺開發(fā)的中國大氣污染物和溫室氣體人為源排放清單模型,主要涵蓋10種大氣污染物和溫室氣體(SO2、NOx、CO、NMVOC、NH3、CO2、PM2.5、PM10、BC 和 OC).根據(jù)相關(guān)日變化曲線[66-67]在 MEIC中加入日變化,盡可能降低排放源對模擬結(jié)果產(chǎn)生的影響,使排放源具有可靠性.生物質(zhì)燃燒排放源來自于 NCAR(FINN)的火災(zāi)清單,具有小時時間分辨率和 1km 水平分辨率[68].垂直粉塵通量采用GOCART粉塵排放方案[69]計算.

    1.4 誤差分析

    為驗證模式模擬結(jié)果的準(zhǔn)確性,本文采用平均偏差(MB),均方根誤差(RMSE)和相關(guān)系數(shù)(R)評估模式模擬結(jié)果.

    2 結(jié)果與討論

    2.1 PM2.5模擬結(jié)果

    由圖2可見,6種方案均較好地模擬出了PM2.5濃度隨時間的變化趨勢,然而 1月1日之后,模擬的風(fēng)速過大疊加增大后的湍流擴散系數(shù)導(dǎo)致ACM2R2和 ACM2R5方案顯著低于觀測值,因此,本文計算了1月1日之前、之后、以及整個時段中的白天、夜間和全天的偏差(表1).由表1可以看出,在1月1日之前,YSU、MYJ、MYNN2和 ACM2均較大地高估了 PM2.5濃度,且夜間偏差大于白天;ACM2R2明顯地修正了該高估現(xiàn)象,而ACM2R5則在夜間轉(zhuǎn)為低估.在 1月 1日之后,YSU、MYJ、MYNN2和ACM2在夜間仍顯著高估PM2.5濃度,而ACM2R2和ACM2R5則轉(zhuǎn)為顯著低估.就整個時段而言,YSU、MYJ和 MYNN2高估 PM2.5濃度,夜間高估更為顯著;ACM2白天輕微低估,夜間顯著高估;而ACM2R2和ACM2R5白天和夜間都為低估,其中ACM2R5低估更多.

    表1 各方案PM2.5濃度模擬值與觀測值的對比統(tǒng)計Table 1 Comparison of simulated and observed PM2.5 concentrations in different schemes

    圖2 各方案PM2.5模擬結(jié)果與觀測結(jié)果對比Fig.2 Comparison of PM2.5 simulation and observation results by thesix schemes

    在整個時段中,YSU、MYJ、MYNN2和ACM2,4種方案對于南京地區(qū) PM2.5濃度模擬結(jié)果與觀測值的MB為24.78~46.04μg/m3,均存在一定的高估,而經(jīng)過修訂后的ACM2R2的MB轉(zhuǎn)為低估為-9.1μg/m3,ACM2R5則進(jìn)一步加大低估,其MB為-17.33μg/m3.YSU、MYJ和MYNN2的均方根誤差約為60μg/m3,而ACM2為49.06μg/m3,ACM2R2和ACM2R5則降低為 40.01,36.01μg/m3.各方案對于南京地區(qū) PM2.5濃度模擬結(jié)果與觀測值的 R均在 0.7左右,其中YSU、MYJ和MYNN2低于0.7,分別為0.65、0.66和0.63,ACM2為0.71,而ACM2R2和ACM2R5則提高為0.77和0.76.

    對整個時段而言,YSU、MYJ、MYNN2和ACM2,4種方案中,ACM2的平均偏差MB和均方根誤差RMSE最小,而R最大,說明其在此次污染事件的模擬中對于3個統(tǒng)計指標(biāo)均表現(xiàn)最優(yōu).而另外3種方案中,YSU的均方根誤差最大,MYJ的平均偏差MB最大,MYNN2的R最小.此外,經(jīng)過進(jìn)一步對湍流擴散系數(shù)的修訂,ACM2R2和ACM2R5在ACM2的基礎(chǔ)上進(jìn)一步提升了這3個統(tǒng)計指標(biāo)與觀測值的接近程度,其中平均誤差由24.7μg/m3優(yōu)化為-9.1,-17.33μg/m3,均方根誤差由 49.06μg/m3減小為40.1,36.1μg/m3,R由0.71提升至0.77和0.76.

    2.2 2m溫度與10m風(fēng)速

    近地層溫度與風(fēng)速是地-氣相互作用中重要的大氣下邊界信息,其體現(xiàn)了地面過程與大氣過程之間相互作用的結(jié)果,其中2m溫度和10m風(fēng)速通常作為代表性指標(biāo)量[70-71].由圖3(a)可以看出,6種方案均很好地模擬出 2m溫度的日變化特征,即夜間溫度低而白天溫度高,午后溫度達(dá)到最大值.模式在白天和夜間對于溫度的模擬均存在不同程度上的低估,白天 6種方案模擬結(jié)果差異較小,偏差在 0.5℃左右;夜間除MYJ方案與觀測數(shù)據(jù)較為一致外其他4種方案偏差均較大,清晨時分偏差可超過2℃,其中在2016年12月31日和2017年1月1日清晨達(dá)將近3℃.

    圖3 2m溫度和10m風(fēng)速的模擬值與觀測值對比Fig.3 Comparison of 2m temperature and 10m wind speedsimulation and observation

    由表 2可知,各方案的平均偏差均為負(fù)數(shù),其中MYJ偏差最小為-0.3℃,而經(jīng)過修訂的ACM2R5偏差最大達(dá)-1.65℃,這是增加的湍流擴散更利用夜間近地層大氣與地表進(jìn)行熱交換,從而使得降溫更為劇烈.均方根誤差中仍是 MYJ最小為 1.48℃,ACM2R5最大為 1.98℃.因為溫度變化趨勢的模擬與觀測十分吻合,6種方案的R都接近1,MYJ為0.97,其余皆為0.96.

    表2 2m溫度與10m風(fēng)速觀測值與模擬值的對比統(tǒng)計Table 2 Statistical comparison results between simulation and observationfor 2m temperature and 10m wind speed

    由圖3(b)可以看出,本次污染過程10m風(fēng)速較小且呈現(xiàn)明顯的日變化特征,白天風(fēng)速較夜間大,中午風(fēng)速最大一般超過1m/s但最大不超過2m/s,夜間風(fēng)速極小,至午夜時分基本接近或低于儀器的相應(yīng)閾值(0.1m/s).6種邊界層方案對10m風(fēng)速的模擬均存在顯著高估,偏差最大可超過2m/s.前人研究也發(fā)現(xiàn)模式傾向于高估城市中的近地層風(fēng)速,其主要原因是城市建筑物對風(fēng)場的減弱作用在模式中被低估[72-73].

    由表2可以看出,各方案10m風(fēng)速的平均偏差均為正值,其中 MYJ偏差最大達(dá) 1.71m/s,ACM2R5偏差最小為 1.11m/s.均方根誤差也是 MYJ最大為1.82m/s,ACM2R5最小為 1.26.R則是 MYJ最小為0.32,ACM2R5最大為 0.37.所以這3個統(tǒng)計量對于2m溫度和10m風(fēng)速剛好相反,6種方案中MYJ模擬得到的溫度偏高風(fēng)速偏大,而ACM2R5則溫度偏低風(fēng)速偏小.這主要是由于各方案得到的邊界層內(nèi)湍流擴散系數(shù)大小所導(dǎo)致.6種方案中 MYJ的湍流擴散系數(shù)最小,不利于大氣向夜間較冷地面散熱,從而溫度稍高;另外較小的湍流擴散系數(shù)也不利于動量向下混合,從而風(fēng)速較高;具有較大湍流擴散系數(shù)的ACM2R5的機制則剛好相反,從而模擬得出較低溫度和較低風(fēng)速.

    2.3 垂直廓線數(shù)據(jù)對比驗證

    污染物濃度最高值出現(xiàn)在12月31日中午,為檢驗不同邊界層方案在垂直方向上的模擬狀況,選取12月30日以及12月31日當(dāng)?shù)貢r間08:00和20:00的2km高度內(nèi)的探空廓線數(shù)據(jù)與6種方案的模擬值進(jìn)行對比,包括溫度、緯向風(fēng)(U風(fēng)),經(jīng)向風(fēng)(V風(fēng))和水汽混合比.

    整體來看,6種方案在垂直方向上對溫度的模擬差異不大,主要的差別集中在地面附近.由圖 4可知,南京地區(qū)30~31日08:00和20:00,4個時次近地面附近均存在逆溫,而6種邊界層方案均較好地模擬出了溫度隨高度的變化趨勢,但由于垂直分辨率較粗,所以并不能捕捉其中的溫度層結(jié)細(xì)節(jié),特別是早上8:00,400m左右高度處和晚上20:00,800m左右高度處的溫度廓線.這里是殘余層與自由大氣的交界處,亦受到平流影響,存在氣象要素的突變.

    圖4 溫度垂直廓線Fig.4 Temperature vertical profile

    對比ACM2和ACM2R2及ACM2R5可知,將模式底部的湍流擴散系數(shù)設(shè)定一個較小的閾值(5,2m2/s),僅對溫度廓線模擬值產(chǎn)生細(xì)微影響,主要的區(qū)別仍集中在模式底層,ACM2R5的 2m溫度的統(tǒng)計量較ACM2偏離觀測值更多.

    由圖 5可以看出,模式對水汽混合比廓線的模擬雖然與觀測值較接近,但存著高估或低估現(xiàn)象.相比于風(fēng)速和溫度,不同邊界層方案之間濕度模擬值的差別較大,最大可達(dá)約 1倍.如圖 5(b)中 30日20:00,900m 高度處,MYNN2、ACM2和 ACM2R5的水汽混合比模擬值是YSU和MYJ方案的近1倍,而圖 5(d)中 31日 20時 600m 處 YSU、ACM2和ACM2R5是MYNN2方案的將近1倍.這些高估低估的轉(zhuǎn)變以及不同方案之間的巨大差別說明除邊界層方案外的其他物理過程(如陸面過程)對水汽的模擬亦有較大影響[74].

    圖5 水汽混合比垂直廓線Fig.5 Vertical profile of water vapor mixing ratio

    圖 6中,分探空的風(fēng)速觀測值垂直分辨率約為300m.與溫度廓線類似,6種邊界層方案的差別仍較小.在地面附近,U風(fēng)和V風(fēng)的風(fēng)速模擬值都傾向于大于觀測值(圖6c除外),但隨著高度的增加,伴隨著風(fēng)向的轉(zhuǎn)變,風(fēng)速模擬值有可能開始小于觀測值(如 400~1200m,1000m高度以上).總體上看,6種方案均較好地模擬出了U風(fēng)和V風(fēng)的廓線變化狀況,只是風(fēng)向轉(zhuǎn)折的高度和風(fēng)速大小與觀測值之間仍存在些許差別.

    圖6 U風(fēng)和V風(fēng)垂直廓線Fig.6 U Wind and V Wind vertical profile

    2.4 邊界層高度對比

    邊界層高度是模式模擬污染物垂直分布的重要參數(shù),其隨時間變化而變化.不同邊界層參數(shù)化方案診斷邊界層高度的方法不盡相同.YSU方案所計算的邊界層高度診斷為達(dá)到臨界理查森數(shù)時所對應(yīng)的高度(穩(wěn)定情況下臨界理查森數(shù)為 0.25,不穩(wěn)定情況下為 0),MYJ方案則將邊界層高度診斷為湍流動能強度下降到臨界值 0.01m2/s2的高度,MYNN2方案的邊界層高度由湍流動能和虛位溫廓線共同決定,ACM2方案與YSU類似采用臨界理查遜數(shù)方法確定邊界層高度(臨界理查森數(shù)為0.25).

    基于激光雷達(dá)探測大氣邊界層高度的方法主要有4種:分別為梯度法、標(biāo)準(zhǔn)偏差法、曲線擬合法和小波協(xié)方差變換法.由激光雷達(dá)接收到的米散射回波信號變化可得激光雷達(dá)距離平方校正信號P(z)z2.梯度法是直接根據(jù) P(z)z2隨高度衰減速度率大小作為判斷大氣邊界層高度的依據(jù).標(biāo)準(zhǔn)偏差法是利用一定高度范圍內(nèi)P(z)z2測量數(shù)據(jù)的離散程度來確定大氣邊界層高度.曲線擬合法是通過構(gòu)建曲線函數(shù)來擬合目標(biāo)曲線,選取最為吻合的曲線所對應(yīng)的高度.小波協(xié)方差變換法則通過探測信號階躍變化的高度來確定大氣邊界層高度[75-76].Fan等[42]利用激光雷達(dá)采用立方根梯度算法計算了這一期間的邊界層高度變化,并利用秒探空所測得邊界層高度對其進(jìn)行了驗證.在夜間,由于殘余層的存在,白天的污染物仍可能懸停半空,造成激光雷達(dá)所測大氣邊界層高度可能存在一定程度的高估.

    將 Fan等[42]利用激光雷達(dá)回波梯度法得到的邊界層高度與模式模擬結(jié)果對比(圖7).由圖7可知,邊界層高度呈現(xiàn)明顯的日變化,一般中午~午后最高,可達(dá)1.5km,夜間最低,可低至200m.期間20日夜間和29日早晨由于殘余層頂部污染物濃度梯度較大,激光雷達(dá)測得的邊界層高度出現(xiàn)跳躍式變化,極大的高估邊界層高度,這主要是由于夜間殘余層頂部被激光雷達(dá)錯誤地識別為邊界層高度所導(dǎo)致.整體而言,6種方案均較好地模擬出了邊界層高度逐日變化以及日變化特征.如12月28日和29日的觀測和模擬的日最大邊界層高度均較高(超過 1km),而12月31日、1月1日和3日則較低(低于800m).另外,從日出時分邊界層開始發(fā)展至中午邊界層發(fā)展到最高的時段,邊界層高度發(fā)展速率的模擬在 6種方案中都較為準(zhǔn)確.然而,傍晚時分邊界層高度下降在模式中通常十分突然(MYNN2除外),這主要是由于地面附近逆溫層開始發(fā)展,而邊界層高度的計算則由逆溫層產(chǎn)生前的混合層頂部高度(約 1km)轉(zhuǎn)變?yōu)榈孛娓浇鏈貙拥母叨?約 100~300m).MYNN2方案中邊界層高度的計算合并了虛位溫廓線的影響,使得傍晚時分邊界層高度的下降較為平緩,與激光雷達(dá)觀測更為接近.夜間,相比于激光雷達(dá)的邊界層高度,模式輸出的邊界層高度均偏低,甚至經(jīng)常低至幾十米.模式結(jié)果中這些夜間過低的邊界層高度模擬值與夜間偏高的污染物濃度模擬值高度相關(guān),但是由于邊界層高度在模式中只是一個診斷量,它同樣只是模式的一個模擬結(jié)果,并不能作為污染物濃度高估的原因.例如,MYNN2傍晚時分邊界層高度下降速率明顯低于其他方案,但它的污染物濃度上升速度卻與其他方案相同;再如夜間通常MYJ和MYNN2的邊界層高度較高,但同樣是這種方案的污染物濃度較大.需進(jìn)一步挖掘不同邊界層方案 PM2.5濃度模擬偏差的原因.整體上看 6種方案在模擬邊界高度方面偏低;其中MYNN2方案模擬最好,平均偏差MB和均方根誤差RMSE最小,R最大,為0.76;其他4種邊界層方案差別不大,R均在0.61左右.

    表3 模擬的邊界層高度誤差統(tǒng)計Table 3 Statistical table of errors between simulated and actual boundary layer heights

    圖7 不同參數(shù)化方案邊界層高度模擬結(jié)果與測量結(jié)果對比Fig.7 Comparison of simulation and measurement results of boundary layer height with different parameterization schemes

    2.5 湍流擴散系數(shù)

    在WRF-Chem模式中,湍流引起的熱量和物質(zhì)擴散作用被視作相同,即熱量湍流擴散系數(shù)和物質(zhì)湍流擴散系數(shù)相等,簡稱為湍流擴散系數(shù).湍流擴散系數(shù)是直接影響污染物擴散的重要參量,它表征了污染物在垂直方向上混合的強度.湍流擴散系數(shù)越大,混合越強,垂直方向上的污染物濃度越均勻,從而地表的污染物濃度越低,反之亦然.

    圖8為6種方案模擬所得的湍流擴散系數(shù)從模式底層到 1.3km 高度內(nèi)的變化.從逐日變化的角度上看,白天湍流擴散系數(shù)較大值(大于 5m2/s)所能達(dá)到的高度與圖7中的日邊界層最大高度很好地相互對應(yīng),白天湍流擴散系數(shù)較大值能達(dá)到 1km左右的高度,12月28日和29日較高,12月31日、1月1日和 3日則較低.從日變化的角度上看,湍流擴散系數(shù)在夜間各高度上都較小,一般不超過 10m2/s;清晨日出后大約 9:00開始快速增大且大值區(qū)快速向上延伸,于 13:00~14:00達(dá)到最大值,最大超過 100m2s-1,最高達(dá) 1km左右.各邊界層方案模擬的湍流擴散系數(shù)在白天和晚上均有較大差異,白天YSU較大值向上延伸的高度最高,但 MYJ和 MYNN的值最大,ACM2、ACM2R2和ACM2R5的值最小.夜間各方案的差異更為明顯,4種未經(jīng)修改的方案在這一時段的湍流擴散系數(shù)常為最低閾值(小于 0.1m2/s),即湍流被強烈抑制,擴散作用微弱.其中,非局地方案YSU和 ACM2在清晨和傍晚均能模擬出湍流擴散系數(shù).YSU傾向于將擴散模擬在100~500m的高度,這主要是由YSU中邊界層頂?shù)膴A卷項所貢獻(xiàn).ACM2則將擴散模擬在貼地 100~200m,這與其模式參數(shù)化方程有關(guān),ACM2的湍流擴散系數(shù)并非像YSU一樣直接使用 K理論,而是將模式各層之間的交換進(jìn)行參數(shù)化,其算法與地面熱通量高度關(guān)聯(lián),在地面通量為正值的情況下其交換系數(shù)通常不會低至閾值.局地方案MYJ和MYNN在清晨時分能夠模擬出較小的交換系數(shù),但在傍晚時分各層交換系數(shù)均觸及模式最小閾值.

    圖8 各方案模擬的低空湍流擴散系數(shù)分布對比Fig.8 Comparison of low altitudeturbulent diffusion coefficient distribution simulated bythe schemes

    由圖9可見,湍流擴散系數(shù)在9:00左右開始升高,于13:00~ 14:00達(dá)最大值,19:00之后則接近0值,直至次日8:00.其中,MYJ和MYNN2整層平均值較大,YSU次之,ACM2最小.

    結(jié)合污染物濃度模擬(圖 10)可以發(fā)現(xiàn),各方案白天湍流擴散系數(shù)差別雖然較大,但是污染物濃度模擬卻無顯著差別;而晚上各方案污染物濃度有所差別,主要是體現(xiàn)在ACM2和另外3種方案之間.傍晚時分 ACM2方案中邊界層依然存在一定的擴散,使得污染物濃度上升較慢,從而午夜時分其污染物濃度值顯著低于另外3種方案;而YSU方案中傍晚時分雖然也模擬出擴散,但它位于夾卷層,對于地面污染物的擴散并不起作用.由圖 9可以看出,對于PM2.5濃度的垂直分布,清晨時分(8:00)ACM2與其他方案并無顯著區(qū)別,但再傍晚時分(20:00)ACM2的地面污染物濃度顯著低于另外3種方案.

    圖9 邊界層高度內(nèi)各方案的平均湍流擴散系數(shù)Fig.9 Averaged turbulent diffusion coefficient within atmospheric boundary layer for the schemes

    對于ACM2R2和ACM2R5方案,在1月1日之前,它們的 PM2.5濃度模擬與觀測較為接近,由于加強了夜間的湍流擴散,模式夜間高估PM2.5濃度的現(xiàn)象得到顯著改善,但是1月1日的天氣形勢造成當(dāng)日地面通量較小,各方案湍流擴散系數(shù)皆較小(圖9),且各方案在當(dāng)日清晨極大地低估了 PM2.5濃度下降的速率(圖 2),在這種情況下,2或 5m2/s的下限閾值似乎已經(jīng)過大,使得當(dāng)日 ACM2R2和 ACM2R5中PM2.5濃度低至20 μg/m3,并于后期皆處于低估PM2.5濃度的狀態(tài).由于2,5m2/s湍流擴散系數(shù)的設(shè)置僅僅是示例湍流擴散系數(shù)對污染物濃度的影響,而并非提出一個新的閾值設(shè)定,所以未進(jìn)一步對閾值的大小進(jìn)行調(diào)試.湍流擴散系數(shù)的大小應(yīng)主要依據(jù)觀測進(jìn)行調(diào)整,而并非僅僅依據(jù)模擬結(jié)果調(diào)整閾值,需將針對各邊界層方案的特點,依據(jù)觀測數(shù)據(jù)在參數(shù)化過程中調(diào)整傍晚及夜間時分的湍流擴散系數(shù)算法,從而達(dá)到優(yōu)化污染物夜間擴散模擬的目的.另外,2,5m2/s湍流擴散系數(shù)閾值對氣象要素的影響極小,這主要是因為氣象要素同時受到地面和邊界層上層自由大氣背景環(huán)流的共同影響.而PM2.5則主要來自地面排放,PM2.5濃度的垂直差異顯著,擴散對其影響巨大.

    圖11 各方案的湍流擴散系數(shù)高度Fig.11 Height variation of turbulent diffusion coefficientfrom the different schemes

    2.6 討論

    湍流擴散系數(shù)的不同可能是6種方案PM2.5濃度模擬差異的重要原因.本文中 6種方案模擬的湍流擴散系數(shù)在白天和晚上均有較大差異,其中在夜間的差異對PM2.5濃度模擬的影響顯著.非局地方案YSU和ACM2在傍晚均能模擬出低層的湍流擴散系數(shù),但YSU傾向于將擴散模擬在100~500m的高度,這主要是由YSU中邊界層頂?shù)膴A卷項所貢獻(xiàn);而ACM2則將擴散模擬在貼地100~200m之內(nèi),這更利于源于地表的污染物向上擴散,所以 ACM2方案在夜間相比其他方案模擬所得的地表PM2.5濃度更低.此外,提高模式低層湍流擴散系數(shù)最小閾值的ACM2R2和ACM2R5方案則在模擬時段前3d基本修正了ACM2的PM2.5濃度高估,但在后期由于擴散過強產(chǎn)生低估.后續(xù)研究中我們將針對各邊界層方案的特點,修正傍晚及夜間時分湍流擴散系數(shù)算法,使得污染物夜間擴散得模擬更為合理.

    Hu等[77]使用WRF-Chem模擬臭氧時發(fā)現(xiàn)湍流擴散系數(shù)的大小可能是影響夜間臭氧濃度模擬的重要因子.王繼康等[23]認(rèn)為,較好的垂直擴散條件會使臭氧濃度模擬增強.董春卿等[27]指出,不同邊界層參數(shù)化方案和湍流輸送的模擬差異,可能是影響近地面PM2.5模擬差異的主要原因.上述研究均表明湍流擴散系數(shù)對污染物濃度的模擬存在顯著影響.但除此之外模式模擬所得的 PM2.5濃度和多種因子有關(guān),如本文采用的所有邊界層方案均模擬得到較低的邊界層高度和溫度,較低的溫度在水汽混合比相當(dāng)?shù)那闆r下可導(dǎo)致較高的相對濕度,從而對PM2.5的化學(xué)過程產(chǎn)生影響.此外,排放源的準(zhǔn)確性、化學(xué)過程刻畫的偏差、其他物理過程參數(shù)化的偏差、以及再分析資料氣象場的誤差等都是模擬誤差的可能來源.在后續(xù)研究中,將需既有針對性的逐一探討,也需綜合地考慮它們的共同作用,從而使得模擬結(jié)果進(jìn)一步接近觀測.

    3 結(jié)論

    3.1 6種方案均較好地模擬出了PM2.5濃度隨時間的變化趨勢,其中YSU、MYJ、MYNN2和ACM2,4種邊界層方案在夜間對 PM2.5濃度的模擬均存在不同程度上的高估,而修訂了湍流擴散系數(shù)閾值的ACM2R2和ACM2R5則顯著地降低高估,甚至轉(zhuǎn)為低估.從平均偏差,均方根誤差和相關(guān)系數(shù) 3個統(tǒng)計量上看,ACM2R2和ACM2R5最接近觀測值,ACM2次之,而YSU的均方根誤差最大,MYJ的平均偏差最大,MYNN2的相關(guān)系數(shù)最小.

    3.2 6種方案模式對2m溫度、10m風(fēng)速和風(fēng)溫濕廓線的模擬與觀測均有所偏差,但模擬結(jié)果均顯示了較為合理的日變化特征及隨高度變化的趨勢.不同邊界層方案之間氣象要素差異并不是 PM2.5濃度模擬差異的主要原因.

    3.3 各方案均較好地模擬出了邊界層高度逐日變化以及日變化特征,但是在傍晚時分邊界層高度下降在模式中通常時分突然(MYNN2除外),這是由于邊界層高度的計算由逆溫層產(chǎn)生前的混合層頂部高度(約1km)轉(zhuǎn)變?yōu)榈孛娓浇鏈貙拥母叨?約100至300m)所導(dǎo)致.MYNN2方案中邊界層高度的計算合并了虛位溫廓線的影響,使得傍晚時分邊界層高度的下降較為平緩.夜間,6種邊界層方案得到的邊界層高度均偏低,甚至經(jīng)常低至幾十米.模式結(jié)果中這些夜間過低的邊界層高度模擬值與夜間偏高的污染物濃度模擬值高度相關(guān).

    猜你喜歡
    邊界層擴散系數(shù)湍流
    基于HIFiRE-2超燃發(fā)動機內(nèi)流道的激波邊界層干擾分析
    重氣瞬時泄漏擴散的湍流模型驗證
    基于Sauer-Freise 方法的Co- Mn 體系fcc 相互擴散系數(shù)的研究
    上海金屬(2015年5期)2015-11-29 01:13:59
    FCC Ni-Cu 及Ni-Mn 合金互擴散系數(shù)測定
    上海金屬(2015年6期)2015-11-29 01:09:09
    一類具有邊界層性質(zhì)的二次奇攝動邊值問題
    非時齊擴散模型中擴散系數(shù)的局部估計
    非特征邊界的MHD方程的邊界層
    “青春期”湍流中的智慧引渡(三)
    “青春期”湍流中的智慧引渡(二)
    弱分層湍流輸運特性的統(tǒng)計分析
    国产视频首页在线观看| 亚洲av中文字字幕乱码综合| 午夜亚洲福利在线播放| 只有这里有精品99| 欧美一级a爱片免费观看看| 在线观看美女被高潮喷水网站| 汤姆久久久久久久影院中文字幕| 99热这里只有是精品50| 在线免费十八禁| 校园人妻丝袜中文字幕| 大码成人一级视频| 精品人妻视频免费看| 亚洲国产高清在线一区二区三| 久久亚洲国产成人精品v| 99热6这里只有精品| 啦啦啦中文免费视频观看日本| 美女国产视频在线观看| 男女边摸边吃奶| 老司机影院毛片| 看黄色毛片网站| 狠狠精品人妻久久久久久综合| 国产精品久久久久久精品电影小说 | 22中文网久久字幕| 在线播放无遮挡| 日韩视频在线欧美| 国产在线一区二区三区精| 免费看av在线观看网站| 精品一区二区三区视频在线| 最近中文字幕高清免费大全6| 大码成人一级视频| 免费av毛片视频| 免费观看a级毛片全部| 男女边摸边吃奶| 22中文网久久字幕| 国产免费福利视频在线观看| 国产精品久久久久久久电影| 久久久久久久久久久免费av| 久久久久久久久久成人| 国产成人午夜福利电影在线观看| 精品久久久噜噜| 高清在线视频一区二区三区| 嘟嘟电影网在线观看| 久久久久九九精品影院| 国语对白做爰xxxⅹ性视频网站| 久久影院123| 国产一区二区亚洲精品在线观看| 99久久中文字幕三级久久日本| av国产免费在线观看| 纵有疾风起免费观看全集完整版| 国产黄片美女视频| 久久精品国产亚洲av涩爱| 国产精品人妻久久久久久| 免费人成在线观看视频色| 国产精品麻豆人妻色哟哟久久| 欧美成人午夜免费资源| 国产精品一二三区在线看| 国产黄色视频一区二区在线观看| 日韩 亚洲 欧美在线| 91精品伊人久久大香线蕉| 99九九线精品视频在线观看视频| 精品一区在线观看国产| 九色成人免费人妻av| 日韩三级伦理在线观看| 色综合色国产| 可以在线观看毛片的网站| 日韩欧美一区视频在线观看 | 亚洲欧洲国产日韩| 免费看a级黄色片| 2022亚洲国产成人精品| 男女啪啪激烈高潮av片| 国产永久视频网站| 欧美精品一区二区大全| 午夜福利网站1000一区二区三区| 国产熟女欧美一区二区| 老司机影院毛片| 成年女人看的毛片在线观看| 熟女电影av网| 欧美日韩一区二区视频在线观看视频在线 | 最近最新中文字幕大全电影3| 亚洲美女搞黄在线观看| 91精品国产九色| 18禁动态无遮挡网站| 亚洲精品国产av蜜桃| 嘟嘟电影网在线观看| av网站免费在线观看视频| 中文字幕久久专区| 97在线人人人人妻| 中文字幕人妻熟人妻熟丝袜美| 涩涩av久久男人的天堂| 免费少妇av软件| 国产色婷婷99| 日韩 亚洲 欧美在线| 18禁裸乳无遮挡免费网站照片| 中国三级夫妇交换| 久久久成人免费电影| 少妇丰满av| 亚洲av成人精品一二三区| 亚洲精品成人久久久久久| 国产一区二区在线观看日韩| 久久久久精品久久久久真实原创| 国产成人91sexporn| 亚洲国产精品专区欧美| 精品人妻视频免费看| av免费观看日本| 国内揄拍国产精品人妻在线| 日本wwww免费看| 亚洲精品乱久久久久久| 亚洲av成人精品一二三区| 久久精品久久久久久噜噜老黄| 国产精品秋霞免费鲁丝片| 丰满人妻一区二区三区视频av| 久久久精品94久久精品| 亚洲欧洲日产国产| 国产精品三级大全| 深爱激情五月婷婷| 国产黄片视频在线免费观看| 亚洲aⅴ乱码一区二区在线播放| 亚洲欧美日韩卡通动漫| 69人妻影院| 性色avwww在线观看| 午夜爱爱视频在线播放| 久久久久网色| 在线免费观看不下载黄p国产| 久久精品久久久久久久性| 国产亚洲一区二区精品| 超碰97精品在线观看| 免费看日本二区| 亚洲久久久久久中文字幕| av在线app专区| 成年女人在线观看亚洲视频 | 精华霜和精华液先用哪个| 日韩精品有码人妻一区| 久久久久精品久久久久真实原创| 爱豆传媒免费全集在线观看| 国产一区亚洲一区在线观看| 色哟哟·www| 久久久精品欧美日韩精品| 久久97久久精品| 精品久久久精品久久久| 高清欧美精品videossex| 丝袜脚勾引网站| 成年免费大片在线观看| 极品教师在线视频| 在线亚洲精品国产二区图片欧美 | 国产在线一区二区三区精| 一级爰片在线观看| 亚洲在久久综合| 五月开心婷婷网| 成人亚洲欧美一区二区av| 身体一侧抽搐| 日韩av在线免费看完整版不卡| 国产精品偷伦视频观看了| 久久99热这里只有精品18| 亚洲一级一片aⅴ在线观看| 男女边吃奶边做爰视频| 18禁在线播放成人免费| 一区二区三区四区激情视频| 天天躁日日操中文字幕| 亚洲国产精品成人久久小说| 久热久热在线精品观看| 丝袜喷水一区| 欧美变态另类bdsm刘玥| 亚洲最大成人中文| 国产精品久久久久久久电影| 国产亚洲5aaaaa淫片| 街头女战士在线观看网站| 久久久久久久久大av| 亚洲最大成人av| 欧美成人一区二区免费高清观看| 亚洲欧美日韩另类电影网站 | 国产免费视频播放在线视频| 亚洲精品日本国产第一区| 如何舔出高潮| 久久久久久久亚洲中文字幕| 熟女av电影| 一个人看的www免费观看视频| 两个人的视频大全免费| 欧美日韩精品成人综合77777| 亚洲av免费在线观看| 视频中文字幕在线观看| 午夜福利视频1000在线观看| 美女视频免费永久观看网站| 久久久久久久亚洲中文字幕| 国产毛片在线视频| 亚洲av日韩在线播放| 婷婷色综合大香蕉| 中文天堂在线官网| 全区人妻精品视频| av网站免费在线观看视频| 26uuu在线亚洲综合色| 有码 亚洲区| 午夜福利视频精品| 下体分泌物呈黄色| 欧美极品一区二区三区四区| av在线老鸭窝| 搞女人的毛片| 在线看a的网站| 欧美日韩精品成人综合77777| 久久久久久国产a免费观看| 久久久久久久久大av| 欧美性猛交╳xxx乱大交人| 国产成人精品一,二区| 极品少妇高潮喷水抽搐| 男女边吃奶边做爰视频| 日本猛色少妇xxxxx猛交久久| 国产视频首页在线观看| 欧美xxxx黑人xx丫x性爽| 亚洲成人久久爱视频| 成人亚洲精品av一区二区| 国产色婷婷99| 国产探花极品一区二区| 亚洲av国产av综合av卡| 九九久久精品国产亚洲av麻豆| xxx大片免费视频| 亚洲国产精品成人综合色| 在线免费十八禁| 国产 一区 欧美 日韩| 亚洲,欧美,日韩| 91久久精品电影网| 亚洲国产精品专区欧美| 亚洲最大成人av| av国产久精品久网站免费入址| 国产精品99久久99久久久不卡 | 七月丁香在线播放| 国产精品精品国产色婷婷| 日韩免费高清中文字幕av| 婷婷色av中文字幕| 亚洲天堂国产精品一区在线| 在线播放无遮挡| 欧美区成人在线视频| 婷婷色综合大香蕉| 综合色av麻豆| 国内精品宾馆在线| 国产69精品久久久久777片| 亚洲图色成人| 久久这里有精品视频免费| 女人久久www免费人成看片| 亚洲精品自拍成人| 久久99热这里只有精品18| 欧美日韩视频高清一区二区三区二| 97热精品久久久久久| 亚洲国产色片| 欧美成人午夜免费资源| 精品久久久久久久人妻蜜臀av| 免费av不卡在线播放| 午夜福利在线在线| 色吧在线观看| 久久ye,这里只有精品| 六月丁香七月| 中文乱码字字幕精品一区二区三区| 女的被弄到高潮叫床怎么办| 另类亚洲欧美激情| 99久久九九国产精品国产免费| 亚洲av中文av极速乱| 亚洲av免费在线观看| 欧美成人一区二区免费高清观看| 天堂俺去俺来也www色官网| 联通29元200g的流量卡| 大陆偷拍与自拍| 2021天堂中文幕一二区在线观| 99热国产这里只有精品6| 偷拍熟女少妇极品色| 久久久色成人| 丝瓜视频免费看黄片| 日日摸夜夜添夜夜添av毛片| 国产av码专区亚洲av| 91午夜精品亚洲一区二区三区| 日韩精品有码人妻一区| 国产一区二区亚洲精品在线观看| 最近手机中文字幕大全| 欧美成人精品欧美一级黄| 精品亚洲乱码少妇综合久久| 99re6热这里在线精品视频| 纵有疾风起免费观看全集完整版| 少妇被粗大猛烈的视频| 又粗又硬又长又爽又黄的视频| 69人妻影院| 麻豆乱淫一区二区| 丝袜喷水一区| 色综合色国产| 亚洲av中文av极速乱| 国产一区亚洲一区在线观看| 免费观看性生交大片5| 亚洲精品日本国产第一区| 欧美 日韩 精品 国产| 国产在线男女| 日本-黄色视频高清免费观看| 日韩不卡一区二区三区视频在线| 久久人人爽av亚洲精品天堂 | 综合色丁香网| 99热全是精品| 另类亚洲欧美激情| 边亲边吃奶的免费视频| 肉色欧美久久久久久久蜜桃 | 亚洲精品国产av成人精品| 嘟嘟电影网在线观看| 免费av观看视频| 1000部很黄的大片| 中国三级夫妇交换| av卡一久久| 蜜桃久久精品国产亚洲av| 国产毛片在线视频| 亚洲一区二区三区欧美精品 | 中文字幕免费在线视频6| 一级毛片aaaaaa免费看小| 天堂中文最新版在线下载 | 嫩草影院新地址| 超碰97精品在线观看| 午夜免费男女啪啪视频观看| 1000部很黄的大片| 一边亲一边摸免费视频| 插阴视频在线观看视频| 内地一区二区视频在线| 人妻一区二区av| 久久女婷五月综合色啪小说 | 久久久亚洲精品成人影院| 最近最新中文字幕免费大全7| 一级片'在线观看视频| 久久久久九九精品影院| 伊人久久国产一区二区| 亚洲最大成人中文| 街头女战士在线观看网站| 观看美女的网站| 99热这里只有是精品在线观看| 精品视频人人做人人爽| 纵有疾风起免费观看全集完整版| 国产精品国产三级国产av玫瑰| 性色avwww在线观看| 极品少妇高潮喷水抽搐| 精品国产一区二区三区久久久樱花 | 国内精品宾馆在线| 国产乱人偷精品视频| 99视频精品全部免费 在线| 日本wwww免费看| 一区二区三区乱码不卡18| 亚洲精品一区蜜桃| 九九在线视频观看精品| 美女国产视频在线观看| 国产在视频线精品| 免费少妇av软件| 麻豆久久精品国产亚洲av| xxx大片免费视频| 国产精品久久久久久精品电影小说 | 少妇猛男粗大的猛烈进出视频 | 国产免费福利视频在线观看| 国产精品人妻久久久影院| 春色校园在线视频观看| 在线亚洲精品国产二区图片欧美 | 精品国产一区二区三区久久久樱花 | 国产在线一区二区三区精| 成人亚洲欧美一区二区av| 黄片wwwwww| 性插视频无遮挡在线免费观看| 九九爱精品视频在线观看| 人人妻人人澡人人爽人人夜夜| 成人国产av品久久久| 美女cb高潮喷水在线观看| 99热这里只有是精品50| av在线蜜桃| 欧美少妇被猛烈插入视频| 国产精品国产三级国产专区5o| 久久精品国产亚洲av涩爱| 亚洲国产av新网站| 99久久精品国产国产毛片| 久久久欧美国产精品| eeuss影院久久| 欧美激情在线99| 观看免费一级毛片| 午夜福利高清视频| 久久人人爽av亚洲精品天堂 | 日韩一区二区三区影片| 国产午夜精品一二区理论片| 日本与韩国留学比较| 久久韩国三级中文字幕| 国产综合懂色| 欧美精品一区二区大全| av在线天堂中文字幕| 欧美极品一区二区三区四区| 国产高清不卡午夜福利| 亚洲精品色激情综合| 少妇熟女欧美另类| 精品酒店卫生间| 日韩制服骚丝袜av| 亚洲va在线va天堂va国产| 免费人成在线观看视频色| 欧美亚洲 丝袜 人妻 在线| 亚洲色图综合在线观看| 新久久久久国产一级毛片| 国产永久视频网站| 少妇裸体淫交视频免费看高清| 久久久久久久国产电影| 亚洲成人av在线免费| 麻豆国产97在线/欧美| 一级av片app| 久久久久久久久久人人人人人人| 高清日韩中文字幕在线| 国产 精品1| 欧美+日韩+精品| 久久久久久伊人网av| 乱系列少妇在线播放| 欧美人与善性xxx| 国产成人免费观看mmmm| 亚洲人与动物交配视频| 久久鲁丝午夜福利片| 丰满人妻一区二区三区视频av| 午夜爱爱视频在线播放| 久久鲁丝午夜福利片| 一本色道久久久久久精品综合| 久久人人爽人人片av| 搞女人的毛片| 少妇人妻 视频| 国产淫片久久久久久久久| 高清视频免费观看一区二区| 国产精品久久久久久精品电影| 久久久久九九精品影院| 国产精品人妻久久久久久| 日韩制服骚丝袜av| 国产成人aa在线观看| 黄色配什么色好看| 九色成人免费人妻av| 晚上一个人看的免费电影| 波野结衣二区三区在线| 国产亚洲午夜精品一区二区久久 | 边亲边吃奶的免费视频| 欧美日韩视频高清一区二区三区二| 亚洲久久久久久中文字幕| 身体一侧抽搐| 日本av手机在线免费观看| 国产精品麻豆人妻色哟哟久久| 汤姆久久久久久久影院中文字幕| 欧美成人一区二区免费高清观看| 久久久久久久亚洲中文字幕| 免费高清在线观看视频在线观看| 一二三四中文在线观看免费高清| 成人综合一区亚洲| 亚洲av欧美aⅴ国产| 久久6这里有精品| 午夜免费观看性视频| 又大又黄又爽视频免费| 国产精品一区www在线观看| 国产精品一区二区三区四区免费观看| 啦啦啦在线观看免费高清www| 天堂俺去俺来也www色官网| 国产av码专区亚洲av| a级毛片免费高清观看在线播放| 超碰97精品在线观看| 欧美日韩一区二区视频在线观看视频在线 | 免费看av在线观看网站| 69人妻影院| 国产 一区 欧美 日韩| 2021天堂中文幕一二区在线观| 亚洲精品国产av成人精品| 国产在线男女| 一区二区av电影网| 国产成人aa在线观看| 看免费成人av毛片| 一级片'在线观看视频| 狠狠精品人妻久久久久久综合| 欧美精品一区二区大全| 免费观看无遮挡的男女| 天堂网av新在线| 欧美日韩视频精品一区| 99久国产av精品国产电影| 国产探花极品一区二区| 国产成人福利小说| 国产久久久一区二区三区| av在线观看视频网站免费| 久久久久九九精品影院| 91狼人影院| 亚洲精品乱久久久久久| 免费看av在线观看网站| 肉色欧美久久久久久久蜜桃 | 一级爰片在线观看| 国产真实伦视频高清在线观看| 国产老妇伦熟女老妇高清| 亚洲av二区三区四区| 国国产精品蜜臀av免费| 九九爱精品视频在线观看| 91午夜精品亚洲一区二区三区| av卡一久久| 高清视频免费观看一区二区| 精品午夜福利在线看| 18禁动态无遮挡网站| 亚洲精品视频女| 禁无遮挡网站| 亚洲最大成人中文| 亚洲精品国产成人久久av| 久久精品国产亚洲av天美| 亚洲激情五月婷婷啪啪| 少妇猛男粗大的猛烈进出视频 | 日韩不卡一区二区三区视频在线| 麻豆成人午夜福利视频| 国产成人精品福利久久| 国产精品三级大全| 亚洲图色成人| 男人舔奶头视频| 18禁裸乳无遮挡免费网站照片| 免费看光身美女| 精品久久国产蜜桃| 久久精品国产亚洲av天美| 国产极品天堂在线| 禁无遮挡网站| 国产成人freesex在线| 日本猛色少妇xxxxx猛交久久| 内地一区二区视频在线| av播播在线观看一区| 最近的中文字幕免费完整| 大话2 男鬼变身卡| 天天躁日日操中文字幕| 特级一级黄色大片| 九色成人免费人妻av| 久久久久久久久久久丰满| 免费观看a级毛片全部| 国产黄色视频一区二区在线观看| 全区人妻精品视频| 国产免费视频播放在线视频| 亚洲无线观看免费| 精品久久久久久久久av| 中文在线观看免费www的网站| 日韩精品有码人妻一区| 亚洲精品456在线播放app| 蜜臀久久99精品久久宅男| 国产精品一区www在线观看| 亚洲国产高清在线一区二区三| 国产成人免费无遮挡视频| 在线精品无人区一区二区三 | 国产精品三级大全| 国产日韩欧美在线精品| 免费av毛片视频| 精品亚洲乱码少妇综合久久| 国产精品av视频在线免费观看| 成年免费大片在线观看| 亚洲精品,欧美精品| 少妇人妻精品综合一区二区| 交换朋友夫妻互换小说| 男插女下体视频免费在线播放| 欧美区成人在线视频| 好男人视频免费观看在线| 国产精品.久久久| 91在线精品国自产拍蜜月| 欧美激情国产日韩精品一区| 色视频在线一区二区三区| 人妻一区二区av| 深爱激情五月婷婷| 亚洲精品一二三| av在线播放精品| 国产探花极品一区二区| 欧美日韩视频高清一区二区三区二| 80岁老熟妇乱子伦牲交| 久久久久久久午夜电影| 日韩电影二区| 男女边吃奶边做爰视频| 国产欧美日韩精品一区二区| 2018国产大陆天天弄谢| 精品久久久噜噜| 日韩三级伦理在线观看| 国产 一区 欧美 日韩| 中国三级夫妇交换| 老女人水多毛片| 亚州av有码| 欧美97在线视频| av又黄又爽大尺度在线免费看| 极品少妇高潮喷水抽搐| 99热网站在线观看| 免费黄频网站在线观看国产| 免费av毛片视频| 精品人妻一区二区三区麻豆| av免费在线看不卡| 观看免费一级毛片| 久久久久久久久大av| 国产精品麻豆人妻色哟哟久久| 成人国产av品久久久| 特大巨黑吊av在线直播| 神马国产精品三级电影在线观看| 男人爽女人下面视频在线观看| 亚洲第一区二区三区不卡| 免费观看的影片在线观看| 日本免费在线观看一区| 一本一本综合久久| 国产精品一区二区性色av| 国产在视频线精品| 制服丝袜香蕉在线| 久久久久精品性色| 我的老师免费观看完整版| 亚洲欧美一区二区三区国产| 在线观看av片永久免费下载| 99热这里只有是精品在线观看| 亚洲成人中文字幕在线播放| 久久精品夜色国产| 亚洲精品久久久久久婷婷小说| 韩国av在线不卡| 日日啪夜夜撸| 国产视频内射| 亚洲成人中文字幕在线播放| 中文字幕免费在线视频6| 亚洲精品成人久久久久久| 色视频www国产| 婷婷色综合大香蕉| 久久久久久久久久久丰满| 亚洲av男天堂| 精品久久久噜噜| 国产亚洲91精品色在线| 国产探花在线观看一区二区| 日本猛色少妇xxxxx猛交久久| 天堂网av新在线| 国产成人91sexporn| 国产成人a区在线观看| 国产成人精品一,二区| 亚洲综合色惰| 观看美女的网站| 免费av观看视频| 亚洲av电影在线观看一区二区三区 | 男人和女人高潮做爰伦理| 日韩av不卡免费在线播放| 日韩一区二区三区影片|