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

邊界層對流對示蹤物抬升和傳輸影響的大渦模擬研究

2015-12-14 09:14:40王蓉黃倩田文壽張強張健愷桑文軍
大氣科學 2015年4期
關鍵詞:效率

王蓉 黃倩 田文壽 張強 張健愷 桑文軍

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

2 甘肅省人工影響天氣辦公室,蘭州730020

3 中國氣象局干旱氣象研究所甘肅省干旱氣候變化與減災重點實驗室,蘭州730020

1 引言

沙塵暴是一種危害極大的災害性天氣,會引起一系列的生態(tài)與環(huán)境問題,如荒漠化、土壤肥力下降、空氣污染等,對人類生命和財產(chǎn)安全都造成了嚴重的危害(趙思雄和孫建華,2013)。沙塵氣溶膠的遠距離輸送,如中國及中亞地區(qū)的沙塵,在合適的條件下可以長距離傳輸?shù)綎|亞的韓國、日本(Tan et al., 2012;Kim et al., 2013),使沙塵暴成為區(qū)域性和全球性的環(huán)境問題。另外,沙塵暴在氣候變化中也扮演著重要角色。礦物沙塵氣溶膠作為地球氣候系統(tǒng)的重要組成之一(IPCC, 2007),一方面,沙塵氣溶膠的輻射效應直接影響輻射收支平衡,進而引起區(qū)域乃至全球的氣候變化(Haywood et al., 2005);另一方面,它可以為冰云的形成提供凝結核,從而影響云微物理結構、光學特性及降水形成(Field et al., 2006),對氣候變化產(chǎn)生間接的效應。另外,沙塵氣溶膠還為浮游生物提供了必需的礦物元素(Mahowald et al., 2005),維持了自然界的生物鏈。

國內外學者對沙塵暴的形成機理、遠距離輸送以及其對輻射和氣候影響等方面已經(jīng)做了大量工作,并取得了一定的研究成果(Yang et al., 2013;O’Loingsigh et al., 2014)。作為沙塵輸送基礎的起沙過程也一直是科學工作者研究的重要內容(李耀輝和張書余,2007;趙琳娜等,2007,Li and Zhang,2014),因為起沙的定量和準確描述是模擬和預報沙塵濃度的基礎(張宏生和李曉嵐,2014)。風蝕起沙是起沙機制的核心內容,影響風蝕起沙的因子有天氣和氣候條件、土壤特性、地表特征和實際土地利用(Shao,2008),其中風和大氣邊界層結構是影響風蝕起沙的關鍵因子,而邊界層對流對邊界層結構和發(fā)展又有重要影響(Huang et al., 2010;黃倩等,2014)。近年來邊界層對流對沙塵抬升和垂直傳輸?shù)挠绊懯艿皆絹碓蕉嗟年P注(Huang et al.,2010;Bozlaker et al., 2013;Ramaswamy, 2014)。尤其是在極端干旱的沙漠地區(qū),其夏季晴天邊界層厚度可以發(fā)展到4~6 km(Gamo, 1996;Marsham et al.,2008a),深厚的邊界層對流一方面可以將沙塵傳輸?shù)捷^高的高度(Takemi, 1999),另一方面為沙塵在水平方向的遠距離輸送提供了有利條件(Iwasaka et al., 2003)。Takemi et al.(2006)模擬研究了中緯度沙漠地區(qū)晴天條件下,由邊界層干對流和積云對流引起的沙塵抬升和傳輸?shù)膭恿^程。研究結果表明,邊界層干對流對沙塵在邊界層內的垂直混合有重要作用,而邊界層干對流和積云對流的耦合能使沙塵從邊界層傳輸?shù)阶杂纱髿狻akemi et al.(2005)的研究還指出,極端干旱的沙漠地區(qū)深厚的干對流活動有利于高空水平動量的下傳,使地表風速增加,從而有助于地表沙塵的抬升。Knippertz et al.(2009)認為,撒哈拉西部夜間低空急流造成的動量下傳有利于第二天邊界層對流的發(fā)展和沙塵的抬升。而Todd et al.(2008)對撒哈拉不同測站的觀測研究也證明了邊界層低層風速的增加對于沙塵抬升的重要作用。另外,Huang et al.(2010)還利用大渦模式,診斷分析了非均勻的地表熱通量對撒哈拉沙漠地區(qū)示蹤物抬升效率及垂直傳輸?shù)挠绊懀芯拷Y果表明由于非均勻地表熱通量引起局地地表風速增大,從而增加示蹤物的抬升效率。另外,非均勻地表熱通量引起的局地環(huán)流有利于沙塵從混合層向撒哈拉殘留層的傳輸。雖然這些研究結果加深了我們對邊界層對流對沙塵抬升及垂直傳輸影響的理解,但是目前系統(tǒng)地分析不同形式的邊界層對流對沙塵抬升和垂直傳輸影響的研究較少,另外,也缺少描述地表熱通量和風切變對沙塵抬升效率和傳輸高度影響的定量關系,而這對于深入研究邊界層對流對沙塵抬升效率和垂直傳輸具有重要意義。

敦煌位于我國西北干旱區(qū),氣候干燥,又與我國沙塵天氣的高發(fā)區(qū)南疆盆地接壤,為該地區(qū)沙塵天氣的形成提供了必要的物質條件,是我國河西走廊沙塵天氣的高發(fā)區(qū)。敦煌地區(qū)白天較強的地表加熱能力和夜間冷卻能力造成該地區(qū)超厚對流邊界層的發(fā)展及演變(張強等,2007;Zhang et al., 2011)。趙建華等(2011)用熱力數(shù)值模型對西北干旱地區(qū)對流邊界層高度的定量分析也表明感熱是西北干旱區(qū)深厚對流邊界層形成的主要原因。本文將在這些研究的基礎上,以敦煌地區(qū)加密觀測期間的實測資料為背景,利用大渦模式,分析干旱區(qū)地表感熱通量(以下稱為地表熱通量)和地轉風切變(以下稱為風切變)對邊界層對流強度及對流形式的影響,并進一步研究不同形式的邊界層對流對示蹤物抬升效率和垂直傳輸?shù)挠绊懀詈蠖康亟o出示蹤物抬升效率及傳輸高度隨地表熱通量和風切變變化之間的關系。

2 模式及方法介紹

本文所利用的英國氣象局大渦模式[Large Eddy Model (LEM) Version 2.4 (Gray et al., 2001)]是一個高分辨率、非靜力平衡的數(shù)值模式,可以用來模擬范圍廣泛的湍流尺度和云尺度的問題(對模式的具體描述見黃倩等(2014))。本研究中模式高度取為6 km,水平區(qū)域為10 km×10 km,水平方向采用等距的網(wǎng)格,為200 m,垂直方向采用隨高度變化的網(wǎng)格距,其中最小的格距在近地面,約1.4 m,最大的在3 km之上,約為158 m。模式中采用了周期側邊界條件和鋼性上下邊界條件,為了減少由模式上界引起的重力波反射,在模式高度約4 km(約為模式高度2/3處)以上應用了牛頓阻尼吸收層。模式中使用的地表地轉風是由NCEP–NACR 機構的 2.5°×2.5°再分析資料計算得到的,地轉風切變是用小球探空資料1 km高度的風速和地表地轉風資料求得。

本文使用的觀測資料是“西北干旱區(qū)陸氣相互作用野外觀測實驗”加密觀測期間2000年6月3日敦煌站的位溫、比濕探空廓線,以及敦煌雙墩子戈壁站的地表熱通量觀測資料。大氣穩(wěn)定度是表征湍流發(fā)展的一個重要參數(shù),以湍流能量為基礎的理查遜數(shù)(Richardson number)同時包含了影響層結穩(wěn)定度的熱力和動力因子,其中與熱力因子有關的地表熱通量決定了影響邊界層湍流發(fā)展的浮力的強弱,而與動力因子有關的風切變能把邊界層湍流(對流)組織成不同的形式(Weckwerth et al.,1999)。本文將在標準試驗的基礎上,通過一系列改變地表熱通量強度和風切變大小(具體增加或減小為標準試驗地表熱通量和風切變的倍數(shù)見表1)的敏感性數(shù)值試驗診斷分析不同地表熱通量和風切變條件下的邊界層對流中,示蹤物的抬升效率及傳輸高度。為了研究邊界層對流對物質傳輸?shù)挠绊懀谀J?00 m高度以下加入絕對濃度為100的被動示蹤物。另外,由于研究區(qū)域是極端干旱的沙漠地區(qū),其波恩比值較大,潛熱通量對模擬結果影響不大,因此在敏感性數(shù)值試驗中改變地表熱通量是指改變地表感熱通量的大小。本研究的試驗設計如表1所示,其中C0代表標準試驗(地表熱通量和風切變?yōu)閷崪y結果),表1中的數(shù)字代表各敏感性數(shù)值試驗中地表熱通量和風切變放大(或縮小)為標準試驗中地表熱通量和風切變的倍數(shù)。

表1 敏感性試驗設計。H和S分別代表地表熱通量和風切變,HC0和 SC0分別代表標準試驗地表熱通量和風切變,√代表有試驗,×代表無試驗,C1至C9分別代表改變地表熱通量或風切變的試驗Table 1 The design of the sensitivity tests. H and S are the surface heat flux and wind shear, respectively. HC0 and SC0 represent the surface heat flux and shear from the standard experiment, respectively. Symbol √ indicates the experiment was performed, and × indicates it was not. C1 to C9 are the experiments with different surface heat fluxes or wind shears

3 模擬結果與分析

6月3日的實測資料和模擬結果都顯示[黃倩等(2014)的圖1],07時[北京時間(BT),下同],有200 m厚的貼地逆溫層,其上是厚度約800 m、強度約為0.02°C m-1的覆蓋逆溫層,10時由于地表受太陽輻射加熱形成的熱泡向上發(fā)展,對流邊界層(CBL)的厚度約為300 m,逆溫層之上是一層清晰可辨的厚度約為3 km的近中性分層的殘留層。14時CBL頂已經(jīng)到達1100 m的高度,而且覆蓋逆溫的強度也有所減弱。到 16時邊界層對流運動將混合層與殘留層貫通為一層厚度約為4 km的超厚邊界層。本文的標準試驗是以 07時的探空廓線為初始場,模擬的不同時次的位溫廓線與實測廓線基本一致。

3.1 不同形式的邊界層對流

邊界層對流有不同的形式,如無組織的對流泡、有組織的邊界層對流卷等(Etling and Brown,1993),而地表熱通量和風切變對邊界層對流的強度和結構有重要影響(Tian et al., 2003; Shin and Hong, 2013),因此本文通過改變模式地表熱通量和風切變大小形成不同形式的邊界層對流。圖1是不同敏感性數(shù)值試驗模擬的不同時次的邊界層位溫廓線。從圖1看出,11時到14時,CBL不斷增暖,而且 CBL的厚度也在逐漸增大。當風切變一定,隨著地表熱通量的增大(試驗 C0、C2、C3、C4),由于近地層熱泡獲得的能量增多,上沖的高度也增高,從而使CBL變暖且厚度增大。如14時試驗C4(地表熱通量增大為標準試驗地表熱通量的1.2倍)的CBL的平均位溫約是316 K,而試驗C2(地表熱通量是標準試驗地表熱通量的0.5倍)的CBL平均位溫只有308 K左右。當?shù)乇頍嵬恳欢ǎS著風切變的增大(試驗C0、C6、C7、C8),CBL也增暖變厚,另外,還注意到當?shù)乇頍嵬恳欢ǎ龃箫L切變時CBL的厚度增長較快,如圖1d中試驗C8的CBL厚度約為1.1 km,而試驗C4的CBL厚度僅有0.7 km。這主要是因為增大風切變加強了夾卷層的夾卷效率,使更多自由大氣的暖空氣向下混合從而使CBL增暖,而且增強的夾卷作用有助于混合層與覆蓋逆溫層的貫通,從而使 CBL增厚(黃倩等,2014)。

圖1 不同試驗模擬的(a)11時(北京時間,下同)、(b)12時、(c)13時、(d)14時的平均位溫廓線。單位:KFig. 1 Simulated profiles of averaged potential temperature (K) for different cases at (a) 1100 BT (Beijing time, the same below), (b) 1200 BT, (c) 1300 BT, and (d) 1400 BT

下面進一步分析不同形式邊界層對流的空間分布特征。圖2是不同敏感性數(shù)值試驗模擬的 13時300 m高度的垂直速度水平剖面圖。從圖2可以看出,風切變不變,減小地表熱通量(試驗C2),邊界層對流的強度明顯減弱(300 m高度氣流的最大上升速度只有2.0 m s-1);而增大地表熱通量(試驗C4),邊界層對流的強度增強(300 m高度氣流的最大上升速度可以達到4.5 m s-1)。當?shù)乇頍嵬坎蛔儯瑴p小風切變(試驗C6),邊界層對流的形式明顯改變,而且對流的強度也較標準試驗的略增強;增大風切變,邊界層對流的強度有所減弱,這與風切變影響湍流渦旋的方向,從而改變邊界層對流的強度有關(Paugam et al., 2010)。另外,增大風切變,邊界層對流出現(xiàn)較弱的有組織的對流卷信號,如圖2g上升氣流和下沉氣流排列較規(guī)則,若將風切變和地表熱通量分別改變?yōu)闃藴试囼烇L切變的2.0倍,地表熱通量改變?yōu)闃藴试囼灥乇頍嵬康?0.8倍(如圖2h),邊界層對流的上升和下沉氣流排列更規(guī)則,但是對流的強度較弱(上升氣流的最大速度為 1.6 m s-1),這與 Weckwerth et al.(1999)的研究結果一致,即邊界層對流卷出現(xiàn)在地表熱通量一定(不小于50 W m-2),風切變較大的條件下。另外,圖2b和圖2e分別為只有浮力驅動和只有動力驅動邊界層對流的試驗模擬的垂直速度,可以看出只有浮力和只有動力驅動的邊界層對流泡都比較細碎,但是浮力驅動的邊界層對流強度明顯大于動力驅動的邊界層對流強度。Moeng and Sullivan(1994)的研究結果也顯示:地表熱通量的大小對邊界層對流的強弱有重要影響,而不同大小的風切變會將邊界層對流組織成不同的形式。

圖2 13 時 300 m 高度的垂直速度(m s-1)水平分布。其中,(a)、(b)、(c)、(d)、(e)、(f)分別為試驗 C0、C1、C2、C4、C5、C6、C8、C9 的結果Fig. 2 Horizontal distribution of vertical velocity (m s-1) at 300 m at 1300 BT for experiments (a) C0, (b) C1, (c) C2, (d) C4, (e) C5, (f) C6, (g) C8, and (h) C9

為了進一步分析不同形式邊界層對流在垂直方向的空間結構,圖3給出了圖2中各敏感性數(shù)值試驗模擬的13時垂直速度垂直剖面圖。從圖3可以看出,當風切變不變時,增大地表熱通量的試驗中(圖3d,試驗C4)邊界層對流發(fā)展的高度(約1 km)大于減小地表熱通量試驗中(圖3c,試驗C2)邊界層對流的高度(約 0.5 km)。當?shù)乇頍嵬坎蛔儠r,增大風切變(圖3g,試驗 C8)邊界層對流發(fā)展的高度減小(約 0.5 km),但是上升氣流和下沉氣流排列更規(guī)則。另外,從圖3還可以看出,只有熱力驅動的邊界層對流中,對流泡垂直向上發(fā)展,如圖3b中的熱泡與地表基本垂直,而隨著風切變的加入,對流泡發(fā)生傾斜,如圖3a中y=0 km、圖3g中y=-3.9 km、y=-1.0 km及y=4.5 km處,這一特點隨著風切變的增大更顯著(圖3h),這與黃倩等(2014)的研究結果一致。

圖3 13 時沿 x=0 m 的垂直速度(m s-1)y-z剖面圖,其中(a)、(b)、(c)、(d)、(e)、(f)分別為試驗 C0、C1、C2、C4、C5、C6、C8、C9 的結果Fig. 3 y–z cross sections of vertical velocity (m s-1) at x=0 m at 1300 BT for experiments (a) C0, (b) C1, (c) C2, (d) C4, (e) C5, (f) C6, (g) C8, and (h) C9

圖4給出了敏感性數(shù)值試驗模擬的不同形式邊界層對流 13時示蹤物絕對濃度及風矢量隨高度分布的垂直剖面圖。從圖4可以看出當?shù)乇頍嵬繙p小為標準試驗地表熱通量的0.5倍時(試驗C2,圖4c),示蹤物傳輸?shù)母叨燃s為 0.6 km,而將地表熱通量增大到標準試驗的1.2倍時(試驗C4,圖4d),示蹤物傳輸?shù)母叨燃s為1.7 km。當?shù)乇頍嵬繛闃藴试囼灥牡乇頍嵬浚L切變減小為0.5倍的標準試驗風切變時(試驗C6,圖4f),示蹤物可以傳輸?shù)酱蠹s1.0 km的高度,如果把風切變增大到標準試驗的1.5倍(試驗C8,圖4g),示蹤物傳輸?shù)母叨让黠@增大,可以達到約2.2 km的高度。圖4b、4e、4h分別為只有熱力驅動邊界層對流(試驗C1)、只有風切變驅動邊界層對流(試驗C5)和有較強對流卷的試驗(試驗C9)結果,從圖中可以看出,示蹤物在有較強對流卷信號的試驗中傳輸?shù)母叨茸畲螅s 3.2 km),另外,雖然在只有熱力驅動的邊界層對流中混合層的高度(約0.7 km)大于只有風切變驅動試驗的邊界層對流中的混合層高度(約 0.3 km),但是試驗 C5中示蹤物傳輸?shù)母叨瓤梢赃_到1.4 km,而在試驗C1中示蹤物的最大傳輸高度只有0.9 km。另外還注意到,增大風切變試驗中的夾卷層厚度都比增大地表熱通量試驗中的夾卷層厚度大,如圖4e、4g、4h與圖4d中示蹤物濃度大于0.1的藍色區(qū)域的比較。從圖4的分析知道,增大地表熱通量和增大風切變都能使示蹤物傳輸?shù)母叨仍龃螅龃蟮乇頍嵬恐饕窃鰪娏诉吔鐚訉α鞯膹姸龋簿褪窃龃罅松蠜_熱泡的能量,從而使示蹤物隨著強上升氣流傳輸?shù)捷^高的高度;而增大風切變主要是增強了夾卷作用、增大了夾卷層厚度,有利于混合層和覆蓋逆溫層的混合貫通,從而CBL明顯增厚,使示蹤物的傳輸高度增大,而且增大風切變比增大地表熱通量更有利于示蹤物傳輸?shù)捷^高的高度。在所有不同形式的邊界層對流中,有較強邊界層對流卷信號的試驗中示蹤物傳輸?shù)母叨茸罡摺?/p>

圖4 試驗 (a) C0、(b) C1、(c) C2、(d) C4、(e) C5、(f) C6、(g) C8、(h) C9模擬的13時風矢量和示蹤物絕對濃度的垂直剖面圖Fig. 4 y–z cross sections of wind fields and absolute concentration of tracer at 1300 BTfor (a) case C0, (b) case C1, (c) case C2, (d) case C4, (e) case C5, (f)case C6, (g) case C8, and (h) case C9

3.2 不同形式的邊界層對流對示蹤物抬升效率的影響

雖然邊界層對流對沙塵的垂直傳輸有重要影響(Takemi et al., 2006;Huang et al., 2010),但是邊界層對流對沙塵抬升效率的影響也不容忽視(Marsham et al., 2008b)。Gillette(1978)的研究結果顯示當?shù)乇盹L速大于某一臨界值時,沙塵才能從地表被抬升。全球礦物沙塵的貢獻中有35%是由對流泡和對流渦旋引起的(Koch and Renno, 2005)。Marsham et al.(2011)通過定義沙塵抬升潛力,研究了夏季深對流形成的冷池對西非沙塵抬升的影響,模擬結果顯示氣候模式中對對流過程的參數(shù)化導致近18%的沙塵抬升潛力的減少。本文將定量地研究不同強度和不同形式的邊界層對流對示蹤物抬升及傳輸?shù)挠绊憽J紫炔捎肅akmur et al.(2004)給出的公式:

來研究不同形式的邊界層對流對示蹤物抬升效率的影響。(1)式中F代表示蹤物抬升效率,u是10 m高度處的風速,uT是臨界風速,只有u≥uT時計算的示蹤物抬升效率有效。這里需要說明的是(1)式中示蹤物抬升效率的單位是m3s-3。

為了更準確地分析地表熱通量和風切變對示蹤物抬升效率和傳輸高度的影響,利用表1中除了C1、C5、C7和C9四個試驗以外的126個改變地表熱通量和風切變的敏感性數(shù)值試驗結果做進一步的分析。其中風切變取為標準試驗C0的風切變,地表熱通量在試驗C0地表熱通量的0.1倍至1.5倍(最大值為37 W m-2到555 W m-2)之間變化;地表熱通量取為標準試驗C0的地表熱通量,風切變在試驗C0風切變的0.03倍至1.5倍(0.52×10-3s-1至24.5×10-3s-1)之間變化。

在研究不同形式的邊界層對流對示蹤物抬升效率和傳輸高度影響之前,首先分析了臨界風速的大小對示蹤物抬升效率的影響。研究結果(圖略)顯示臨界風速值越大,示蹤物的抬升效率越低,這與(1)式中的理論分析結果是一致的。另外,改變地表熱通量和風切變的大小對示蹤物總體抬升效率隨臨界風速的變化趨勢影響不大。

圖5是用最小二乘法擬合的13時示蹤物抬升效率隨地表熱通量變化的直線,其中符合擬合直線變化規(guī)律的試驗結果用黑色的十字表示,而不符合擬合直線變化規(guī)律的試驗結果用藍色十字表示。從圖5可以看出,風切變小于0.6倍標準試驗風切變時(黑十字代表的試驗),示蹤物的抬升效率隨地表熱通量的增大而增強,也就是當風切變小于10.5×10-3s-1時地表熱通量越大示蹤物的抬升效率也越大,而風切變大于0.6倍標準試驗風切變的試驗中(藍十字代表的試驗)示蹤物的抬升效率與地表熱通量沒有這樣的變化規(guī)律。另外從圖5還可以看出,風切變越小(如風切變?yōu)闃藴试囼烇L切變的0.03倍),示蹤物的抬升效率隨地表熱通量的增加而增大得越快。這主要是因為地表熱通量越大,邊界層對流的強度越大,越有利于上層動量的下傳,從而使近地層風速增大,導致示蹤物的抬升效率增大;而風切變越小,熱力作用越顯著,垂直方向湍流運動增強,也有助于上層動量的下傳。圖2和圖3的分析結果也顯示,增大風切變使邊界層對流的強度減弱,因此只有風切變較小時,示蹤物的抬升效率才隨地表熱通量的增大而增加。當風切變小于標準試驗風切變的0.6倍時,示蹤物抬升效率隨地表熱通量變化的關系式可以表示為

圖5 用最小二乘法擬合的13時示蹤物平均抬升效率隨地表熱通量變化的直線。橫坐標是地表熱通量(H)放大(或縮小)為標準試驗地表熱通量(HC0)的倍數(shù),縱坐標是抬升效率F。其中,實線代表最小二乘擬合直線,黑色十字代表符合擬合直線的試驗結果,而藍色十字代表不符合擬合直線的試驗結果Fig. 5 Straight lines (solid) of least square fitting to the tracer mean lifting efficiency (F) vs surface heat flux changes at 1300 BT. The x-axis represents H/HC0, black (blue) crosses denote the experimental results are consistent (inconsistent) with the lines

其中,F(xiàn)表示示蹤物的抬升效率,H表示地表熱通量,a和b是隨風切變大小變化的系數(shù)。a、b的取值范圍分別為3.65×10-2至6.02×10-2、4.88×10-4至 65.7×10-4。

為了進一步分析風切變對示蹤物抬升效率的影響,用最小二乘法擬合圖5中126個試驗中風切變與示蹤物抬升效率之間的關系,如圖6。從圖6可以看出,13時地表熱通量小于標準試驗地表熱通量的1.2倍(444 W m-2)時,示蹤物的抬升效率隨著風切變的增大反而減小,但是當?shù)乇頍嵬吭黾訛樵囼濩0中地表熱通量的1.5倍(555 W m-2)時,只有當風切變不太大時(小于標準試驗的1.0倍)才滿足示蹤物抬升效率隨風切變的增大而減小的規(guī)律。另外還注意到,地表熱通量越大,示蹤物抬升效率隨風切變的增大遞減越快。其中當?shù)乇頍嵬啃∮跇藴试囼灥乇頍嵬康?.2倍,示蹤物抬升效率隨風切變的變化規(guī)律可以表示為

其中,F(xiàn)表示示蹤物抬升效率,S代表風切變,a1、b1、c1和d1是隨地表熱通量變化的系數(shù)。a1、b1、c1和d1分別在4.5×10-2至19.1×10-2、2.84至16.1、66.7至134.2、2.61×103至13.7×103之間變化。通過圖2和圖3的分析已經(jīng)知道隨著風切變的增大,邊界層對流強度減弱,也就是上升氣流的速度減小,這導致上層動量下傳減弱,近地層風速減小,從而使示蹤物的抬升效率減小。而圖6中顯示當?shù)乇頍嵬繛闃藴试囼灥乇頍嵬康?.5倍,風切變也增大為標準試驗風切變的1.5倍時,示蹤物的抬升效率卻明顯增大,這也許與超強的地表熱通量引起上層動量下傳有關。

圖6 用最小二乘法擬合13時示蹤物平均抬升效率隨風切變變化的曲線。橫坐標是風切變(S)放大(或縮小)為標準試驗風切變(SC0)的倍數(shù),縱坐標是抬升效率F。其中,實線代表最小二乘擬合曲線,黑色十字代表符合擬合曲線的試驗結果,而藍色十字代表不符合擬合曲線的試驗結果Fig. 6 Curves (solid) of least square fitting to the tracer mean uplifting efficiency (F) vs surface heat flux changes at 1300 BT. The x-axis represents S/SC0,black (blue) crosses denote the experimental results are consistent (inconsistent) with the lines

圖7是不同試驗模擬的13時示蹤物平均抬升效率隨地表熱通量和風切變的變化圖。從圖7可以看出當?shù)乇頍嵬啃∮跇藴试囼灥乇頍嵬康?.25倍(462.5 W m-2)時,示蹤物的抬升效率隨著風切變的增大而減小(圖7中熱通量較小的虛線區(qū)域);風切變小于標準試驗風切變的 0.6倍(10.5×10-3s-1)時,即圖7中風切變較小的點線區(qū)域,示蹤物的抬升效率隨著地表熱通量的增大而增大。

3.3 不同形式的邊界層對流對示蹤物傳輸高度的影響

圖2和圖3的結果表明不同大小的地表熱通量和風切變影響邊界層對流的強度和形式,從而使示蹤物傳輸?shù)母叨纫膊幌嗤榱诉M一步定量化研究不同形式的邊界層對流對示蹤物傳輸高度的影響,以不同大小地表熱通量和風切變的敏感性數(shù)值試驗結果為基礎,利用最小二乘法分別確定地表熱通量和示蹤物的平均傳輸高度以及風切變和示蹤物平均傳輸高度的定量關系。另外需要說明的是,文中示蹤物的平均傳輸高度是用示蹤物的平均濃度小于0.1的最大高度確定的。

圖7 不同地表熱通量和風切變的敏感性數(shù)值試驗模擬的13時示蹤物平均抬升效率。點線代表風切變(S)為標準試驗風切變(SC0)的0.6倍的值,虛線為地表熱通量(H)是標準試驗地表通量(HC0)1.25倍的值Fig. 7 Simulated tracer mean uplifting efficiency vs various surface heat fluxes and wind shears at 1300 BT. Dotted line shows S/SC0=0.6, dashed line represent H/HC0=1.25

圖8 是用最小二乘法擬合的126個敏感性數(shù)值試驗模擬的示蹤物平均傳輸高度隨地表熱通量變化的曲線,曲線用以下三次擬合多項式表示:

其中,h(單位:m)是示蹤物的平均傳輸高度,H(單位:W m-2)是地表熱通量,a2、b2、c2和d2是與風切變大小有關的系數(shù)。從圖8看出,當風切變一定時,示蹤物的平均傳輸高度隨地表熱通量的增加而增大,也就是邊界層對流的強度增大,示蹤物傳輸?shù)母叨纫苍龈撸@主要是因為邊界層湍流將較大的熱量向上層大氣輸送,使熱泡獲得了較多的能量而上沖到較高的高度,示蹤物隨著熱泡上升也被傳輸?shù)捷^高的高度。但是,從圖8可明顯看出示蹤物的平均傳輸高度和地表熱通量的大小之間并不是簡單的線性關系。

同樣用最小二乘擬合的方法進一步分析地表熱通量一定時示蹤物的平均傳輸高度隨風切變的變化規(guī)律,示蹤物的平均傳輸高度可以用下式表示:

其中,h(單位:m)代表示蹤物的平均傳輸高度,S(單位:s-1)是風切變,a3、b3、c3和d3是由不同大小的地表熱通量決定的系數(shù)。不同敏感性數(shù)值試驗模擬的示蹤物平均傳輸高度和擬合曲線如圖9所示,其中示蹤物的平均傳輸高度用黑色十字表示。從圖9可以看出當風切變較小時,示蹤物平均傳輸高度隨風切變增大變化不大,有的略有減小,而當風切變較大時,示蹤物平均傳輸高度隨風切變的增大而明顯增加。如地表熱通量為1.2倍標準試驗地表熱通量、不同風切變的試驗中(H/HC0=1.2的擬合曲線),當風切變小于0.5倍標準試驗中風切變時,示蹤物的平均傳輸高度隨風切變的增大而略有減小;而當風切變大于0.5倍標準試驗中風切變時,示蹤物平均傳輸高度隨風切變的增大而增大,尤其是當風切變大于標準試驗中風切變時,示蹤物的平均傳輸高度隨風切變的增大迅速增加。另外還注意到,增大地表熱通量,示蹤物平均傳輸高度隨風切變的增加而增大更快。由圖1的分析知道,增大風切變使夾卷效率增強,不僅導致 CBL增暖,還增大了CBL的厚度,因此示蹤物也隨著CBL的增厚而被傳輸?shù)捷^高的高度。為了解釋地表熱通量越大,示蹤物的平均傳輸高度隨風切變的增加而增大更快的原因,進一步分析了不同地表熱通量對應的增大風切變試驗模擬的 13時平均位溫廓線,如圖10所示,當?shù)乇頍嵬枯^小時,熱泡上沖的高度較低,混合層以上的覆蓋逆溫層的厚度較大,逆溫強度也較大(圖10中的實線,地表熱通量為試驗C0地表熱通量的0.1倍),雖然此時風切變值也較大,但是混合層中較弱的對流以及覆蓋逆溫層的較強逆溫和較大的厚度都不利于通過夾卷作用將自由大氣的暖空氣向下混合,因此當?shù)乇頍嵬枯^小時,雖然增大風切變加強了夾卷作用,但是CBL的厚度仍然較小。增大地表熱通量,邊界層對流的強度增大,混合層的厚度也明顯增大,覆蓋逆溫層的強度和厚度也都相應減小,此時由風切變引起的夾卷作用更容易將自由大氣的暖空氣夾卷向下,從而使CBL增暖增厚,如圖10中紅色點線和紅色虛線(地表熱通量分別是試驗 C0地表熱通量的 0.7倍和1.5倍)所示。也就是說地表熱通量越大,覆蓋逆溫的強度和覆蓋逆溫層厚度都減小,此時增大風切變使 CBL厚度增長顯著,從而使示蹤物傳輸?shù)母叨纫苍龃螅@一結果與黃倩等(2014)的研究結果一致。

圖8 用最小二乘法擬合13時示蹤物平均傳輸高度隨地表熱通量變化的曲線。橫坐標是熱通量放大(或縮小)為標準試驗熱通量的倍數(shù),縱坐標是傳輸高度h,。其中黑色十字代表不同模擬試驗的結果Fig. 8 Curves of least square fitting to the mean transmission heights (h) of tracer transport and surface heat fluxes at 1300 BT. The x-axis represents S/SC0,black crosses represent the simulation results from different experiments

圖9 用最小二乘法擬合13時示蹤物的平均傳輸高度隨風切變變化的曲線。橫坐標是風切變放大(或縮小)為標準試驗風切變的倍數(shù),縱坐標是傳輸高度h。其中黑色十字代表不同數(shù)值試驗的模擬結果Fig. 9 Least-squares fit of curves to the mean heights of tracer transport and wind shears at 1300 BT. The abscissa represents the increased (or decreased)times of wind shears; the ordinate represents the transmission height h. Black crosses represent the simulated results from various cases

最后綜合分析地表熱通量和風切變對示蹤物平均傳輸高度的影響。從圖11看出,風切變和地表熱通量對示蹤物平均傳輸高度的影響與圖8和圖9的分析結果基本一致,即風切變一定,示蹤物的平均傳輸高度隨地表熱通量的增加而增大;地表熱通量一定,風切變較大時,示蹤物的平均傳輸高度隨風切變的增加也增大(如圖11中風切變較大的虛線區(qū)域)。值得注意的是,地表熱通量的增大對示蹤物平均傳輸高度的增加沒有風切變的限制,而示蹤物的傳輸高度隨風切變的增加而增大是在一定的風切變條件下成立的,而這一風切變的臨界值與地表熱通量的大小密切相關。因為增大地表熱通量有利于垂直方向的湍流發(fā)展,從而使示蹤物傳輸?shù)母叨仍龃螅辉龃箫L切變有利于水平方向的湍流產(chǎn)生,而且加強了邊界層頂?shù)臏u旋運動,增強了夾卷作用(Kim et al., 2003),而這種夾卷作用只有在混合層中湍流運動較強時,才能有效地把夾卷的自由大氣中的暖空氣向下混合,從而使 CBL增厚,示蹤物傳輸?shù)捷^高的高度。

圖10 不同地表熱通量對應的增大風切變試驗模擬的13時平均位溫廓線。其中實線、點線和虛線分別代表地表熱通量為標準試驗地表熱通量0.1倍、0.7倍和1.5倍的模擬結果Fig. 10 Simulated profiles of potential temperature from the experiments with various surface heat fluxes corresponding to increasing wind shear at 1300 BT.Solid lines, dotted lines, and dashed lines show H/HC0=0.1, 0.7, 1.5, respectively

圖11 不同地表熱通量和風切變的敏感性數(shù)值試驗模擬的13時示蹤物平均傳輸高度Fig. 11 Simulated tracer mean uplifting efficiency for various surface heat flux and wind shear esperiments at 1300 BT

4 結論

本文利用極端干旱的敦煌地區(qū) 2000年夏季典型晴天6月3日野外試驗觀測的位溫、比濕以及風速探空廓線作為初始化條件,并用實測的地表熱通量驅動大渦模式LEM,通過一系列改變模式地表熱通量和風切變的敏感性數(shù)值試驗模擬了不同形式的邊界層對流,并分析了幾種典型邊界層對流的空間結構特征,在此基礎上進一步研究了不同形式的邊界層對流對示蹤物抬升效率以及傳輸高度的影響,并分別給出了地表熱通量和風切變對示蹤物抬升效率和傳輸高度影響的定量描述。分析得到以下結論:

(1)當風切變一定時,隨著地表熱通量的增大,CBL變暖且厚度增大,邊界層對流的強度也增強;當?shù)乇頍嵬恳欢〞r,隨著風切變的增大,CBL也增暖變厚,但是邊界層對流的強度有所減弱。增大地表熱通量,示蹤物被增強的上升氣流傳輸?shù)捷^高的高度;而增大風切變增強了夾卷作用,夾卷層厚度增大,使示蹤物的傳輸高度也增大,而且比增大地表熱通量示蹤物傳輸高度增加得更大。在所有不同形式的邊界層對流中,有較強邊界層對流卷信號的試驗中示蹤物傳輸?shù)母叨茸罡摺?/p>

(2)當風切變小于 0.6倍標準試驗的風切變(10.5×10-3s-1)時,增大地表熱通量加強了上層動量下傳,使近地層風速增大,從而使示蹤物的抬升效率也增大,而且風切變越小,示蹤物的抬升效率隨地表熱通量的增加而增大得越快。而當?shù)乇頍嵬啃∮跇藴试囼灥乇頍嵬康?.25倍(462.5 W m-2)時,由于增大風切變減弱了邊界層對流的強度,影響上層動量的下傳,使示蹤物的抬升效率減小。

(3)風切變一定時,示蹤物的傳輸高度隨地表熱通量的增加而增大;地表熱通量一定,只有當風切變大于一定值時,示蹤物的平均傳輸高度才隨風切變的增加而增大,而風切變的臨界值取決于地表熱通量的大小。增大地表熱通量有利于垂直方向的湍流發(fā)展,從而使示蹤物傳輸?shù)母叨仍龃螅辉龃箫L切變有利于增強邊界層頂?shù)膴A卷作用,而這種夾卷作用只有在混合層中湍流運動較強時,才能有效地把夾卷的自由大氣中的暖空氣向下混合,從而使CBL增厚,示蹤物傳輸?shù)捷^高的高度。

(References)

Bozlaker A, Prospero J M, Fraser M P, et al. 2013. Quantifying the contribution of long-range Saharan dust transport on particulate matter concentrations in Houston, texas, using detailed elemental analysis [J].Environ. Sci. Technol., 47 (18): 10179–10187.

Cakmur R V, Miller R L, Torres O. 2004. Incorporating the effect of small-scale circulations upon dust emission in an atmospheric general circulation model [J]. J. Geophys. Res., 109 (D7), doi:10.1029/2003JD004067.

Etling D, Brown R A. 1993. Roll vortices in the planetary boundary layer: A review [J]. Bound.-Layer Meteor., 65 (3): 215–248.

Field P R, M?hler O, Connolly P, et al. 2006. Some ice nucleation characteristics of Asian and Saharan desert dust [J]. Atmospheric Chemistry and Physics, 6 (10): 2991–3006.

Gillette D. 1978. A wind tunnel simulation of the erosion of soil: Effect of soil texture, sandblasting, wind speed, and soil consolidation on dust production [J]. Atmos. Environ. (1967), 12 (8): 1735–1743.

Gamo M. 1996. Thickness of the dry convection and large-scale subsidence above deserts [J]. Bound.-Layer Meteor., 79 (3): 265–278.

Gray M E B, Petch J, Derbyshire S H, et al. 2001. Version 2.3 of the Met.Office large eddy model [R]. Met Office (APR) Turbulence and Diffusion Rep, 276pp.

Haywood J M, Allan R P, Culverwell I, et al. 2005. Can desert dust explain the outgoing longwave radiation anomaly over the Sahara during July 2003? [J]. J. Geophys. Res., 110 (D5), doi:10.1029/2004JD005232.

黃倩, 王蓉, 田文壽, 等. 2014. 風切變對邊界層對流影響的大渦模擬研究 [J]. 氣象學報, 72 (1): 100–115. Huang Qian, Wang Rong, Tian Wenshou, et al. 2014. Study of the impacts of wind shear on boundary layer convection based on the large eddy simulation [J]. Acta Meteorologica Sinica (in Chinese), 72 (1): 100–115.

Huang Q, Marsham J H, Parker D J, et al. 2010. Simulations of the effects of surface heat flux anomalies on stratification, convective growth, and vertical transport within the Saharan boundary layer [J]. J. Geophys. Res.,115 (D5), doi:10.1029/2009JD012689.

IPCC. 2007. Fourth Assessment Report, Working Group I Report: The Physical Science Basis, Intergovernment Panel on Climate Change,download at http://.ipcc.ch/report/ar4 [2015-3-24]

Iwasaka Y, Shibata T, Nagatani T, et al. 2003. Large depolarization ratio of free tropospheric aerosols over the Taklamakan desert revealed by Lidar measurements: Possible diffusion and transport of dust particles [J]. J.Geophys. Res., 108 (D23), doi:10.1029/2002JD003267.

Knippertz P, Ansmann A, Althausen D, et al. 2009. Dust mobilization and transport in the northern Sahara during SAMUM 2006–A meteorological overview [J]. Tellus B, 61 (1): 12–31.

Kim H M, Kay J K, Yang E G, et al. 2013. Statistical adjoint sensitivity distributions of meteorological forecast errors of Asian dust transport events in Korea [J]. Tellus B, 65: 20554.

Kim S W, Park S U, Moeng C H. 2003. Entrainment processes in the convective boundary layer with varying wind shear [J]. Bound.-Layer Meteor., 108 (2): 221–245.

Koch J, Renno N O. 2005. The role of convective plumes and vortices on the global aerosol budget [J]. Geophys. Res. Lett., 32 (18), doi:10.1029/2005GL023420.

Li X L, Zhang H S. 2014. Soil moisture effects on sand saltation and dustemission observed over the Horqin sandy land area in China [J]. J. Meteor.Res., 28 (3): 444–452.

李耀輝, 張書余. 2007. 我國沙塵暴特征及其與干旱關系的研究進展 [J].地球科學進展, 22 (11): 1169–1176. Li Yaohui, Zhang Shuyu. 2007.Review of the research on the relationship between sand-dust storm and arid in China [J]. Advances in Earth Science (in Chinese), 22 (11):1169–1176.

Moeng C H, Sullivan P P. 1994. A comparison of shear- and buoyancy-driven planetary boundary layer flows [J]. J. Atmos. Sci., 51(7): 999–1022.

Mahowald N M, Baker A R, Bergametti G, et al. 2005. Atmospheric global dust cycle and iron inputs to the ocean [J]. Global Biogeochemical Cycles,19 (4), doi:10.1029/2004GB002402.

Marsham J H, Parker D J, Grams C M, et al. 2008a. Observations of mesoscale and boundary-layer scale circulations affecting dust transport and uplift over the Sahara [J]. Atmospheric Chemistry and Physics, 8 (23):6979–6993.

Marsham J H, Parker D J, Grams C M, et al. 2008b. Uplift of Saharan dust south of the intertropical discontinuity [J]. J. Geophys. Res., 113 (D21),doi:10.1029/2008JD009844.

Marsham J H, Knippertz P, Dixon N S, et al. 2011. The importance of the representation of deep convection for modeled dust-generating winds over West Africa during summer [J]. Geophys. Res. Lett., 38 (16),doi:10.1029/2011GL048368.

O’Loingsigh T, McTainsh G H, Tews E K, et al. 2014. The Dust Storm Index (DSI): A method for monitoring broadscale wind erosion using meteorological records [J]. Aeolian Research, 12: 29–40.

Paugam R, Paoli R, Cariolle D. 2010. Influence of vortex dynamics and atmospheric turbulence on the early evolution of a contrail [J].Atmospheric Chemistry and Physics, 10 (8): 3933–3952.

Ramaswamy V. 2014. Influence of tropical storms in the northern Indian Ocean on dust entrainment and long-range transport [M]// Tang D L, Sui G J. Typhoon Impact and Crisis Management. Berlin Heidelberg:Springer, 149–174.

Shin H H, Hong S Y. 2013. Analysis of resolved and parameterized vertical transports in convective boundary layers at gray-zone resolutions [J]. J.Atmos. Sci., 70 (10): 3248–3261.

Shao Y P. 2008. Physics and Modelling of Wind Erosion [M]. Dordrecht:Kluwer Academic Publishing, 467pp.

Takemi T. 1999. Structure and evolution of a severe squall line over the arid region in Northwest China [J]. Mon. Wea. Rev., 127 (6): 1301–1309.

Takemi T, Yasui M, Zhou J, et al. 2005. Modeling study of diurnally varying convective boundary layer and dust transport over desert regions [J].Scientific Online Letters on the Atmosphere, 1: 157–160.

Takemi T, Yasui M, Zhou J, et al. 2006. Role of boundary layer and cumulus convection on dust emission and transport over a midlatitude desert area[J]. J. Geophys. Res., 111 (D11), doi:10.1029/2005JD006666.Tan S C, Shi G Y, Wang H. 2012. Long-range transport of spring dust storms in Inner Mongolia and impact on the China seas [J]. Atmos.Environ., 46: 299–308.

Tian W, Parker D J, Kilburn C A D. 2003. Observations and numerical simulation of atmospheric cellular convection over mesoscale topography[J]. Mon. Wea. Rev., 131 (1): 222–235.

Todd M C, Bou Karam D, Cavazos C, et al. 2008. Quantifying uncertainty in estimates of mineral dust flux: An intercomparison of model performance over the Bodélé Depression, northern Chad [J]. J. Geophys.Res., 113 (D24), doi:10.1029/2008JD010476.

Weckwerth T, Horst T, Wilson J. 1999. An observational study of the evolution of horizontal convective rolls [J]. Mon. Wea. Rev., 127 (9):2160–2179.

Yang Y Q, Wang J Z, Niu T, et al. 2013. The variability of spring sand-dust storm frequency in Northeast Asia from 1980 to 2011 [J]. Acta Meteorologica Sinica, 27 (1): 119–127.

張宏升, 李曉嵐. 2014. 沙塵天氣過程起沙特征的觀測試驗和參數(shù)化研究進展 [J]. 氣象學報, 72 (5): 987–1000. Zhang Hongsheng, Li Xiaolan.2014. Review of the field measurements and parameterization for dust emission during sand-dust events [J]. Acta Meteorologica Sinica (in Chinese), 72 (5): 987–1000.

張強, 趙映東, 王勝, 等. 2007. 極端干旱荒漠區(qū)典型晴天大氣熱力邊界層結構分析 [J]. 地球科學進展, 22 (11): 1150–1159. Zhang Qiang,Zhao Yingdong, Wang Sheng, et al. 2007. A study on atmospheric thermal boundary layer structure in extremely arid desert and Gobi region on the clear day in summer [J]. Advances in Earth Science (in Chinese), 22 (11):1150–1159.

Zhang Qiang, Zhang Jie, Qiao Juan, et al. 2011. Relationship of atmospheric boundary layer depth with thermodynamic processes at the land surface in arid regions of China [J]. Science China Earth Sciences, 54 (10):1586–1594.

趙琳娜, 孫建華, 王超, 等. 2007. 2006年春季最強沙塵暴過程的數(shù)值分析 [J]. 氣候與環(huán)境研究, 12 (3): 309–319. Zhao Linna, Sun Jianhua,Wang Chao, et al. 2007. The numerical characteristics of a severe dust storm over North China in the spring of 2006 [J]. Climatic and Environmental Research (in Chinese), 12 (3): 309–319.

趙建華, 張強, 王勝. 2011. 西北干旱區(qū)對流邊界層發(fā)展的熱力機制模擬研究 [J]. 氣象學報, 69 (6): 1029–1037. Zhao Jianhua, Zhang Qiang,Wang Sheng. 2011. A simulative study of the thermal mechanism for development of the convective boundary layer in the arid zone of Northwest China [J]. Acta Meteorologica Sinica (in Chinese), 69 (6):1029–1037.

趙思雄, 孫建華. 2013. 近年來災害天氣機理和預測研究的進展 [J]. 大氣科學, 37 (2): 297–312. Zhao Sixiong, Sun Jianhua. 2013. Study on mechanism and prediction of disastrous weathers during recent years[J]. Chinese Journal of Atmospheric Sciences (in Chinese), 37 (2): 297–312.

猜你喜歡
效率
你在咖啡館學習會更有創(chuàng)意和效率嗎?
提升朗讀教學效率的幾點思考
甘肅教育(2020年14期)2020-09-11 07:57:42
注意實驗拓展,提高復習效率
效率的價值
商周刊(2017年9期)2017-08-22 02:57:49
引入“倒逼機制”提高治霾效率
質量與效率的爭論
跟蹤導練(一)2
提高食品行業(yè)清潔操作的效率
OptiMOSTM 300V提高硬開關應用的效率,支持新型設計
“錢”、“事”脫節(jié)效率低
主站蜘蛛池模板: 久久a级片| 热久久综合这里只有精品电影| 91日本在线观看亚洲精品| 毛片网站免费在线观看| 国产白浆在线| 六月婷婷精品视频在线观看| 国产精品久久久久久久久kt| 久久综合色播五月男人的天堂| 99久久精品无码专区免费| 88av在线看| 在线观看免费黄色网址| 日韩毛片在线视频| 在线精品视频成人网| 亚洲高清资源| 乱人伦中文视频在线观看免费| 精品人妻一区无码视频| 国产男人的天堂| 日韩欧美中文在线| 97超爽成人免费视频在线播放| 日韩中文字幕亚洲无线码| 久久综合婷婷| а∨天堂一区中文字幕| 国产人成在线观看| 国产精品19p| 亚洲aⅴ天堂| 美女黄网十八禁免费看| 香蕉在线视频网站| 欧类av怡春院| 国产色婷婷视频在线观看| 国产欧美日韩另类精彩视频| 亚洲无码高清一区| 国产天天射| 天天色综网| 亚洲日本在线免费观看| 国产永久免费视频m3u8| 亚洲va精品中文字幕| 欧美中文一区| 99久久精品免费看国产电影| 欧美成在线视频| 国产无码在线调教| 日韩性网站| 国产91小视频在线观看| 色婷婷电影网| 一级黄色片网| 日日拍夜夜操| 日韩A级毛片一区二区三区| 国产九九精品视频| 亚洲六月丁香六月婷婷蜜芽| 极品国产一区二区三区| 91无码网站| 亚洲天堂成人| 久久国产成人精品国产成人亚洲 | 日本人妻一区二区三区不卡影院 | 色综合久久88色综合天天提莫 | AV无码无在线观看免费| 国产麻豆精品久久一二三| 久久黄色小视频| 亚洲成人高清无码| 鲁鲁鲁爽爽爽在线视频观看| 久久黄色毛片| 色老二精品视频在线观看| 亚洲乱码精品久久久久..| 毛片免费视频| 国产永久无码观看在线| 中国黄色一级视频| 精品人妻AV区| 久久永久视频| 久久一色本道亚洲| 伊人色天堂| 色婷婷成人网| 久久精品一卡日本电影| 一级片免费网站| 色综合天天操| 天堂中文在线资源| 99热这里只有成人精品国产| 九九久久99精品| 一本大道东京热无码av| 亚洲美女一区二区三区| 人妻精品全国免费视频| 99视频在线免费观看| 成人国产一区二区三区| 亚洲天堂精品视频|