典型非球形冰晶粒子的凝華增長數值模擬試驗

葛 森 孫繼明 牛生杰

?

典型非球形冰晶粒子的凝華增長數值模擬試驗

葛 森1, 2孫繼明1, 2牛生杰1

1南京信息工程大學大氣物理學院/氣象災害預報預警與評估協同創新中心，南京210044,2中國科學院大氣物理研究所云降水物理與強風暴實驗室，北京100029

在凝華增長過程中，冰晶的形狀隨著溫度和濕度的改變而改變，準確模擬冰晶粒子的演變對于提高云模式的模擬能力起著非常重要的作用。在現有的云模式中，冰晶形狀通常假設為球形，而在實際大氣中，冰晶形狀十分復雜。本研究中，我們根據冰晶凝華增長理論模型建立了一個單個冰晶粒子增長模型，模擬了溫度分別為 －1°C～－30°C時，單個典型非球形冰晶粒子的凝華增長過程。與風洞觀測數據相對比，該模型能夠抓住單個冰晶粒子的軸長，質量以及縱橫比隨溫度和濕度的變化過程。我們進一步將該理論增長模型應用到群粒子的凝華增長過程的模擬。我們釆用歐拉二維正定平流輸送法（MPDATA）模擬了典型非球形冰晶群粒子的凝華增長，并對比分析了在不同縱橫比分辨率下的模擬效果以及溫度變化對冰晶形狀的影響，結果表明運用該數值方法可以合理地模擬出群粒子在凝華增長過程中縱橫比的演變。與目前采用的拉格朗日—歐拉混合平流算法比較，該算法能夠耦合到歐拉動力框架下的分檔云模式中去，這對我們研究冰晶粒子形狀對云微物理過程和動力過程的影響，以及它們對冰粒子凝華增長的反饋作用具有非常重要的科學意義。

冰晶 凝華增長 縱橫比 多維正定平流輸送法（MPDATA）

1 引言

冰晶粒子是構成云體結構的主要組分，在混合云降水過程中起著非常重要的作用，并對雷暴閃電的形成以及對大氣輻射平衡的影響等方面有著重要作用（李娟和毛節泰，2006；周廣強等，2006）。與液相粒子相比，冰晶粒子具有復雜多變的形狀，如片狀、柱狀、輻枝狀、柱帽狀、空心柱狀等（Pruppacher and Klett，1997；Bailey and Hallett，2002，2004），并且不同的冰晶形狀會產生不同的下落末速度（Locatelli and Hobbs，1974；Fukuta，1980）、不同的碰撞效率（黃庚等，2007）、不同的光學散射性質（Stephens et al.，1990）、不同的融化效率（龔乃虎，2007），從而對于云的生命史產生不同的影響（Harrington et al., 1999；Morrison et al., 2005；Fridlind et al., 2007；Avramov and Harrington，2010），這些復雜性增加了準確模擬云微物理變化過程的難度?；诒ЯＷ釉谠莆⑽锢磉^程中的關鍵性作用，準確模擬冰晶粒子在云微物理過程中的演化過程對于提高云模式模擬效果具有重要意義。

在冰晶粒子增長演化過程中，凝華增長是極其重要的微物理過程之一。由于采樣困難等原因，冰晶粒子凝華增長過程的研究大多通過設計不同條件的風洞實驗來實現：Fukuta（1969）進行了在－3～－20°C溫度范圍下，小型冰晶粒子凝華增長1min后的形狀特征；Ryan et al.（1974）報道了在－5～－9°C溫度范圍下，增長3 min后，冰晶的形狀特征；Ryan et al.（1976）將其在1974年的冰晶增長研究溫度擴展至－3～－21°C范圍內；Takahashi and Fukuta（1988）（簡稱：TF88）以及Takahashi et al.（1991）（簡稱：TE91）利用改進后的過冷云風洞報道了在－3～－23°C溫度范圍下，增長至30 min時的冰晶增長速度；Fukuta and Takahashi（1999）進一步歸納整理了風洞的實驗結果；Castellano et al.（2014）在云室中觀測了－6°C、－10°C以及－20°C三個溫度下，冰晶增長66～300 s后的形狀變化。伴隨著云室觀測資料的增加，云物理學家也試圖定量化研究冰晶增長的規律[具體的理論總結可見Fukuta and Takahashi（1999）]。目前，大部分模式仍采用球形代替冰晶形狀（Lin et al., 1983；Reisner et al., 1998；Thompson et al., 2004），但由于球形與真實的冰晶形狀相差較遠，會導致涉及冰晶形狀的微物理過程模擬精度降低。因此，一些學者開始考慮將冰晶的非球形形狀加入云模式（Mitchell et al., 1990；Chen and Lamb, 1994b；Harrington et al., 1995；Meyers et al., 1997；Hashino and Tripoli, 2007），但將冰晶的形狀加入歐拉動力框架下的云模式仍是一大難點。

冰晶形狀雖然復雜，但云中冰晶的基本晶型為六角板狀冰晶，一般用軸與軸分別表示冰晶主晶面與棱面的半軸長。在考慮冰晶形狀的模式中，關于晶面和棱面軸長變化的解決方案一般有以下兩種：第一種解決方案為通過冰晶的質量—尺度冪次關系來進行冰晶形狀的模擬（Mitchell et al., 1990；Harrington et al., 1995；Meyers et al., 1997；Woods et al., 2007；Thompson et al., 2008；Morrison and Grabowski，2008，2010）。不同溫度區域內所對應的冰晶形狀具有不同的系數，這導致溫度區域與區域間有間斷性，并且由于該系數為固定觀測環境下得出，數據使用范圍具有很強的局限性，因此在云模式中不具有普遍適用性。第二種解決方案為Chen and Lamb（1994a）基于電容理論提出的“質量分布假設”解決冰晶增長過程中縱橫比的變化，模擬單個冰晶粒子凝華增長參數化的理論模型框架（簡稱：CL94）。基于該理論框架，并將冰晶形狀因子有效地耦合在云物理模型中，進行了多個冰晶粒子凝華增長的數值試驗（Chen, 1992；Chen and Lamb, 1994b；Hashino and Tripoli，2007，2008，2011；Sulia and Harrington，2011；Harrington et al., 2013a）。Harrington et al.（2013b）進行了兩種方案的數值模擬對比，第二種方案更能抓住冰晶的形狀演變的趨勢，并且模擬效果明顯優于第一種方案。然而，在計算冰晶凝化增長的分檔數值平流時，這些數值模擬試驗均采用Chen and Lamb（1994b）提出的拉格朗日—歐拉混合分檔數值平流算法，在拉格朗日框架下，追蹤模擬每個氣塊內冰晶的凝華增長。該算法卻無法耦合到歐拉動力框架下的云模式中。

在本研究中，我們基于CL94模型理論框架，發展建立了一個單個冰晶粒子增長模型，并模擬了單個冰晶粒子在不同溫度下，冰晶形狀在凝華增長過程中的演變。我們進一步將該模型應用到了冰晶群粒子的凝華增長數值試驗，然而我們選用了歐拉二維正定平流輸送算法（MPDATA），從而使我們能夠在歐拉動力框架下，模擬研究與冰晶形狀有關的更為復雜的云微物理過程和動力過程。

2 試驗模式介紹

2.1 單個冰晶粒子凝華增長模式

單個冰晶粒子凝華增長模式基于Chen and Lamb（1994a）（CL94）提出的冰晶粒子凝華增長理論，該模式用于模擬固定環境條件下單個冰晶粒子的凝華增長過程，即在拉格朗日框架下，追蹤單個冰晶粒子在不同的溫度環境下，各個增長參數的演變過程。該理論主要模擬片狀和柱狀冰晶，將冰晶視為橢球狀。CL94的理論框架中主要由兩個基本公式構成。第一個公式是基于電容模式得到的水汽擴散質量增長公式：

柱狀（>）：，（2）

其中，為冰晶縱橫比，即冰晶的縱軸與橫軸之比，＝/，當＞1時，表示冰晶為柱狀，當＜1時，冰晶為片狀，＝1時，則表示冰晶為球狀冰晶。冰晶的水汽擴散質量增長公式是源于傳統的電容模式計算公式，大多云模式中模擬冰晶粒子的質量凝華增長過程均用傳統的計算公式，但一些學者（Nelson，1994；Sulia and Harrington，2011）指出傳統的計算公式最大的限制就是計算冰晶凝華增長時假定冰晶具有恒定的縱橫比，也就是說冰晶的縱橫比在凝華增長過程中是不變的，即軸與軸的軸長之比不變。但在真實的冰晶凝華增長過程中，由于冰晶粒子晶面與棱面增長速度的不同，隨著時間的增加，冰晶粒子的縱橫比會不斷的變化，因此Chen and Lamb（1994a）首次提出了“質量分布假設”解決冰晶增長過程中縱橫比變化的難題：

其中，d/d是指沿軸與軸方向的線性增長速率，、分別為冰晶主晶面、棱面的質量凝華系數，即為水汽在軸方向的凝華增長效率。在公式（4）中，()被稱作為冰晶的內在增長比，是冰晶主晶面與棱面的凝華系數比值，是一個只與溫度相關的數值，由冰晶粒子不同晶型的表面動力學過程所決定。Chen and Lamb（1994a）利用不同的觀測實驗得到的內在增長比，擬合為一條適用于CL94理論框架的內在增長比曲線。雖然目前已經有一些學者利用不同的方法得出了不同的內在增長比（Sei and Gonda，1989；Hashino and Tripoli；2008），但在本研究中仍采用CL94擬合得到的數據，因為，此數據已經被證明在CL94框架中可以較好地模擬出冰晶凝華增長過程中形狀變化（Chen and Lamb，1994a）。隨著冰晶凝華增長，枝狀、空心部分不斷出現，溫度、水汽密度等因素變化，冰晶密度也會隨之發生變化，因此在模式中需要考慮加入冰晶的凝華密度來反映冰晶密度的變化。在本研究中，凝華密度為冰晶粒子凝華增長過程中質量增長部分的密度，采用的是Chen and Lamb（1994a）所給出的凝華密度公式[見文獻Chen and Lamb（1994a）的公式（42）]，該公式是基于Miller and Young（1979）觀測數據而得出的經驗公式。

對于單個冰晶粒子而言，成功準確模擬其在不同溫度下質量、軸長等變化是將冰晶形狀因子嵌入云模式的研究基礎。通過以上設定及公式建立單個冰晶粒子的凝華增長模式，我們將進行在一定環境條件下單個冰晶粒子的凝華增長試驗。

2.2 群粒子冰晶凝華增長試驗

云微物理過程極其復雜，而冰晶粒子的形狀會對碰撞、淞附、融化等過程產生影響。但由于云微物理過程的復雜性以及冰晶粒子形狀對各個微物理過程均有影響，造成在模式中加入冰晶形狀有很多困難。因此，為將冰晶形狀因子加入到歐拉動力框架下的分檔云模式中，我們先在一個只包含凝華過程的理想試驗模型進行檢驗。該試驗模型取自Sun et al.（2012）的云和氣溶膠一維半分檔模式，我們采用了與原分檔云模式相同的質量分檔，即130個水凝物質量檔，90個氣溶膠質量檔，相對應的等球體半徑范圍分別為8.0×10?3～2.4×104μm和8.0×10?3～2.32×102μm。為描述冰晶粒子的數密度分布，用與質量的自然對數有關的數密度分布函數ice(ln) 來表示。同時，ice(ln)dln則表示在ln與lndln之間的單位體積的冰晶粒子數，dln取ln21/2。在質量維度上，用ice(m) 代表ice(ln) dln，是指水凝物質量檔的檔數，從1到130，即i+1=21/2i。并且，該模式中為了追蹤每個水凝物的氣溶膠，因此冰晶粒子的數密度為ice(ln, lnAP)，其中表示水凝物質量，AP表示氣溶膠質量，原模式還對氣溶膠進行了分類，此處暫不考慮。

為將冰晶粒子的形狀加入試驗模型中，增加表征冰晶形狀的縱橫比分檔，也就是說，試驗模型中冰晶粒子的數密度與冰晶質量、氣溶膠質量以及縱橫比大小有關，其表示變為ice(ln, lnAP, l)，其中表示冰晶粒子的縱橫比維度。與選擇ln作為質量分檔不同，在此選取lg作為縱橫比維度的分檔標準，這樣可以在一定程度上保證模式的精度。同時，ice(lg)dlg表示在lg與lg＋dlg之間的單位體積的冰晶粒子數，但在縱橫比維度上，用ice(i)表示ice(lg)dlg。為了評估縱橫比分檔數對模擬效果的影響，本試驗設置兩組不同的縱橫比分檔方案：方案一為41個縱橫比分檔，其中dlg取lg(9/8×21/2)，即i+1＝9/8×21/2i；方案二為73個縱橫比分檔，dlg取log（11/12×21/2），即i+1=11/12× 21/2i，選擇此取值可使兩個方案的縱橫比數值范圍均在1×10?4～1×104。而對于球狀冰晶（＝1），在這兩種分檔方案中的縱橫比檔數分別為第21個檔和第37個檔，根據此檔進行冰晶粒子形狀的劃分，第1～20（或1～36）檔為片狀冰晶，縱橫比范圍為1×10?4＜＜1，而第22～41（或38～73）檔為柱狀冰晶，縱橫比范圍1＜＜1×104，需說明的是，由于取lg作為分檔標準，因此無法做到第21個檔和第37個檔的縱橫比等于1，但此處已盡量使其接近于1，其數值分別為1.07與1.1，此處認為在一定范圍內不影響冰晶粒子形狀的模擬。

因為冰晶粒子形狀因子的加入，需要考慮模式的數值算法問題。為解決冰晶粒子在質量與縱橫比維度的粒子平流計算問題，Chen and Lamb（1994b）提出了一種拉格朗日—歐拉混合分檔計算方法，在以往運用CL94理論添加冰晶粒子形狀因子的分檔云模式中，均采用此數值算法進行模擬。該算法無法運用到歐拉動力框架下的分檔云模式中去。因此，本研究首次提出釆用與以往不同的數值算法進行冰晶粒子形狀演變的模擬。

云和氣溶膠一維半分檔模式（Sun et al., 2012）的數值算法為一維正定平流輸送法，此算法無法進行冰晶粒子濃度在質量與縱橫比維度的同時平流，因此該試驗模式中數值算法選擇Smolarkiewicz et al.（1990）、Smolarkiewicz（2006）的歐拉二維正定平流輸送法（MPDATA）的非振動解方案，該方案是基于通量修正傳輸方案（FCT），為求解平流方程使用迎風差分方案，并且結合反擴散修正。MPDATA具有近似二階近似，在計算過程中，通過迭代方程一步步縮小隱式耗散。MPDATA的基本平流方程為

其中，和是軸和軸方向的流速，是非負標量場。在理想試驗模型中，和即冰晶質量和縱橫比的變化速率，為冰晶粒子的數濃度。模式中具體平流方程為

其中，ice(ln, lnAP, lg) 是冰晶粒子的數密度，在計算平流過程中需要質量和縱橫比隨時間的變 化速率，即dlnd和dlgd。我們根據Chen and Lamb（1994a）提出的理論中可根據水汽擴散凝華增長得dlnd，但需要對縱橫比變化速率進行推導，根據縱橫比與體積以及體積與質量之間的關系式：

從而可以推出冰晶粒子的縱橫比變化速率與質量變化速率之間的關系式：

其中，G()是內在增長比，dep是凝華密度，是冰晶的整體密度。但是由于在該試驗模式中考慮冰晶的凝華密度有一定難度，因此此處假設冰晶在增長過程中密度不變，即凝華密度大小與冰晶的密度相同，這一缺點將在以后加入分檔云模式時加以改進。

釆用此數值算法不但可以解決模式中冰晶粒子的質量與縱橫比演變的數值平流算法問題，同時使得在歐拉動力框架下，將冰晶形狀因子直接耦合到分檔云模式中，研究有關冰晶粒子形狀因子有關的云微物理過程和動力過程成為了可能。

3 模擬結果

3.1 不同初始半徑冰晶粒子的凝華增長

單個冰晶粒子模型用于模擬冰晶粒子凝華增長過程時，我們只研究環境溫度對于冰晶演變的影響。初始環境場條件為：液面飽和度為1，氣壓為1000 hPa，初始冰晶粒子為球形。該模式時間步長為1 s，通過分別模擬環境場為－1～－30°C下冰晶粒子的凝華增長，研究冰晶的形狀隨著溫度變化。并且已有研究表明初始粒子的半徑對于冰晶增長有一定影響，因此對初始半徑不同的冰晶粒子進行模擬對比，初始半徑分別選擇1 μm、5 μm、10 μm、20 μm，而其它環境場條件相同。

圖1表示不同初始半徑的冰晶粒子，在不同溫度下，凝華增長10 min后的軸與軸軸長。一些學者（Libbrecht，2005；盛裴軒等，2003）通過 總結前人所進行的冰晶凝華增長實驗表明：0～ －4°C、－9～－21°C為片狀冰晶增長區，－4～ －9°C、－21～－30°C為柱狀冰晶增長區。本模式模擬結果與冰晶實驗數據吻合，很好地模擬出冰晶形狀隨溫度的變化規律。并且因為初始半徑不同，圖中模擬結果呈現出一定的規律：初始冰晶粒子半徑越小，冰晶增長速率越快（Sheridan et al., 2009）；初始半徑越小的冰晶粒子，越容易增長為極端形狀，如在片狀冰晶增長區，初始半徑為1 μm的粒子經過10 min后，軸數值均比初始半徑大的模擬結果大，軸數值均比其他初始大的模擬結果小，而在柱狀冰晶增長區，則軸的模擬結果也達到最大，軸的模擬結果反而最?。辉冢?°C與－21°C附近，雖然初始半徑不同，但軸與軸的值大小相同，主要因為此時冰晶近乎球形增長。這與Sulia and Harrington（2011）得的模擬結果一致。

圖1 不同初始半徑的冰晶粒子軸長模擬結果與風動實驗數據對比圖。圖中各條線為冰晶粒子模擬10 min后的冰晶粒子的軸長，初始半徑依次是1 μm（實線）、5 μm（折線）、10 μm（點線）、20 μm（點折線），其中黑色粗線為a軸，黑色細線為c軸。點圖為風洞實驗數據：實心黑點為TF88風洞試驗中增長10 min后a軸數據，實心三角為其c軸數據; 空心黑圈為TE91風洞實驗中增長10 min后a軸數據，空心三角為其c軸數據

與Takahashi and Fukuta（1988）和 Takahashi et al.（1991）的風洞實驗結果相比，單個冰晶粒子模式能夠模擬出不同溫度下冰晶軸長的演變。但在－4～－9°C內冰晶粒子的軸軸長模擬結果偏高，同時軸軸長模擬結果偏低，而－9～－21°C間的模擬結果更佳。并且由于在－23～－30°C中冰晶形狀較為復雜，其晶型對于過飽和度的高低極為敏感，也缺乏風洞數據，因此無法對比。但整體而言，模式模擬的冰晶形狀更偏向極端，也就說較風洞中冰晶形狀而言，模擬結果的柱狀（片狀）的軸（軸）軸長數值更大。

圖2顯示了10 min凝華增長后的單個冰晶粒子的質量模擬結果，單個冰晶質量的峰值出現在－15°C、－6°C、－23°C，這與Chen and Lamb（1994a）、Sulia and Harrington（2011）的模擬結果相同。粒子增長速率最快出現在－15°C下，以1 μm粒子為初始粒子的模擬中，可看到10 min的凝華增長后，冰晶粒子質量達到8.4 μg，這與該溫度下觀測值十分接近。在－9°C以及－21°C下，單個冰晶粒子質量增長較慢，這與冰晶在此溫度下近乎球形冰晶的增長特征有關，因為冰晶的不規則形狀會使冰晶的凝華增長速率加快。與實驗數據相對比，模式能夠模擬出冰晶粒子質量在不同溫度下的變化特征，但是與軸長的模擬結果相同，在－4～－9°C下質量模擬結果偏大，在－9～－21°C下，模擬結果更接近于實驗值。

圖2 不同初始半徑的冰晶粒子質量模擬結果與風洞實驗結果對比圖。圖中模擬結果為以初始半徑依次為1 μm（實線）、5 μm（折線）、10 μm（點線）、20 μm（點折線）的增長10 min后冰晶質量變化圖，點圖為風洞實測數據：實心黑點為TF88冰晶質量，空心黑圈為TE91冰晶質量

縱橫比是判斷冰晶形狀的重要依據，圖3為縱橫比的模擬結果與風洞實驗的對比。單個冰晶粒子在－6°C、－15°C以及－23°C溫度下縱橫比出現極值，這表征在此溫度下冰晶粒子易出現極端形狀。在－15°C時，單個冰晶粒子縱橫比達到最低值，1 μm初始粒子經過10 min的凝華增長后縱橫比變為0.006。－6°C時，縱橫比達到最大值82。模式所模擬的縱橫比隨溫度變化特征與觀測數據整體趨勢相似，抓住了冰晶隨著溫度的降低呈現片狀—柱狀—片狀—柱狀的增長趨勢，但模擬結果比觀測結果更易出現極端形狀，如在柱狀冰晶增長區，模擬結果偏大，而在片狀冰晶增長區，模擬結果偏小。相對比而言，片狀冰晶的模擬效果較好，與風洞數據結果相近，柱狀冰晶的模擬結果偏大。

圖3 不同初始半徑的冰晶粒子縱橫比模擬結果與風洞實驗結果對比圖。圖中模擬結果為以初始半徑依次為1 μm（實線）、5 μm（折線）、10 μm（點線）、20 μm（點折線）的增長10 min后冰晶縱橫比變化圖，點圖為風洞實測數據：空心黑圈為TE91冰晶縱橫比

3.2 群粒子凝華增長對比試驗

大量的觀測以及實驗數據（Ryan et al., 1976；Pruppacher and Klett，1997；Bailey and Hallett，2002，2004）表明－15°C以及－6°C為典型的片狀和柱狀冰晶增長溫度，上述單個冰晶粒子的模擬結果與此一致。因此，在群粒子凝華增長試驗中我們選取這兩個溫度進行對比試驗，而將呈現等速球形增長溫度－9°C作為對比參考溫度。

試驗模式為理想條件下一群冰晶粒子的凝華增長，并且因為只考慮了凝華增長過程，環境中的水汽是唯一提供冰晶粒子進行凝華增長的來源，所以模式中所設環境中初始的冰相過飽和度約為37%。同時在初始條件下，根據云中實際情況，設定的群粒子總數濃度為32個 L?1，圖4為方案一中群粒子凝華增長試驗的初始粒子分布情況，粒子半徑分布范圍為0.07～4.5 μm，相對應的質量范圍為2.05×10?9～3.8×10?4，在縱橫比方向，初始粒子均在第21個縱橫比分檔中，而方案二則是初始冰晶粒子均在第37個檔中，即初始粒子均看為球形冰晶，兩個方案中其他設定均相同。通過在試驗模式中進行10 min的凝華增長模擬，對比10 min后兩種方案的冰晶粒子的質量和縱橫比的大小情況。

圖4 群粒子凝華增長初始粒子濃度分布圖。橫坐標（下）為縱橫比，縱坐標（左）為粒子球形等效半徑，橫坐標（上）為縱橫比的分檔，縱坐標（右）為水凝物的質量分檔

圖5初始溫度為－15°C時凝華增長10 min后群粒子濃度分布圖：（a）縱橫比分檔方案一的模擬結果；（b）縱橫比分檔方案二的模擬結果

－6°C是典型的柱狀冰晶增長溫度，因此初始溫度為－6°C時冰晶凝華增長10 min后（圖6），可以看出冰晶粒子均呈柱狀冰晶，甚至為針狀冰晶。這是因為在過飽和條件下，大量觀測以及數據結果表明：－6°C是冰晶增長為針狀冰晶的典型溫度。該溫度下模擬結果與－15°C時結果相似：兩種不同的縱橫比分辨率會導致模擬結果有所差異。圖6a為方案一的模擬結果：冰晶粒子等球體半徑分布范圍為60～200 μm，縱橫比大約在3～800之間，高值區粒子的縱橫比范圍為10～50。圖6b為方案二的模擬結果：冰晶粒子等球體半徑分布范圍為70～150 μm，縱橫比大約在10～400之間，高值區粒子的縱橫比范圍為20～30。兩種方案相比：方案一模擬的冰晶粒子縱橫比、質量分布范圍更廣，高值區粒子數濃度更高。

圖6 同圖5，但為初始溫度為－6°C

在單個冰晶粒子模式的模擬結果中可看到，與－15°C和－6°C相比，冰晶在－9°C下增長速率較慢。這是因為該溫度為冰晶粒子片狀與柱狀的轉換溫度區域，冰晶粒子會傾向于球形增長，而凝華過程中，冰晶粒子的形狀增長傾向對于冰晶粒子增長的速率會有一定影響。根據所得CL94內在增長比數據，－9°C時內在增長比為1.18，即表明在該溫度下冰晶粒子仍有向柱狀增長的趨勢。圖7為初始溫度為－9°C時，群粒子凝華增長10 min后粒子濃度分布，與圖5、圖6相比，冰晶粒子的質量增長相差不多，但是縱橫比有明顯差異，大部分冰晶粒子近乎球形，方案一中縱橫比范圍為：1～13，方案二的縱橫比范圍為1～10。兩個方案相比，方案一的模擬結果中，群粒子的分布范圍也更廣，高值區的數濃度更高。

圖7 同圖5，但為初始溫度為－9°C

CL94理論中冰晶粒子形狀與內在增長比的大小有關，而內在增長比僅與溫度有關，單個冰晶粒子模式中僅考慮一個冰晶粒子的凝華增長，所以凝華潛熱釋放的影響可忽略不計，在一定溫度下，內在增長比不會發生變化。但是在群粒子凝華增長試驗模型中，一群冰晶粒子凝華增長會釋放大量潛熱，環境溫度升高，內在增長比也會跟隨溫度的變化而變化，這會導致在凝華增長過程中，冰晶縱橫比和質量的變化速率會隨著時間改變。由于初始溫度為－6°C的試驗中溫度變化最為顯著，所以選擇分析該溫度下，不同增長時段的群粒子濃度變化作為個例說明。

在初始溫度為－6°C的條件下，冰晶凝華增長10 min后，環境溫度上升了1.52°C，同時內在增長比改變了六次（圖8）。而在柱狀冰晶增長區內，內在增長比降低會導致冰晶縱橫比的增長速率會變慢，同時軸方向的增長減慢，軸方向的增長加快，但是整體上冰晶仍向柱狀冰晶增長。圖中還表示了兩種分檔方案的溫度以及過飽和度的對比：方案二的溫度上升的較高，同時水汽消耗也較大，但數值上兩種方案差距不大。根據內在增長比隨時間的變化，我們選擇1、3、6、10 min的群粒子濃度分布疊加表示其對冰晶縱橫比以及質量變化速率的影響，分檔方案選擇方案二。

圖8 在－6°C下，增長10 min內溫度、過飽和度（冰相）、內在增長比隨時間變化曲線圖。實線和折線為方案一和方案二的模擬結果，黑色細線為溫度變化，黑色粗線為冰相過飽和度變化，黑色實點為方案二內在增長比變化

圖9 初始溫度為－6°C時，增長1、3、6、10 min群粒子濃度分布疊加圖

4 討論

目前，大部分云模式中仍釆用球形冰晶或固定的冰晶形狀，這完全忽略了冰晶粒子形狀對各個微物理過程的影響。將非球形冰晶加入云模式中不僅可以提高模式的模擬精度，而且對于未來提高天氣預報的準確度以及對云微物理過程的研究具有重要意義。在本研究中，我們分別進行了單個冰晶粒子與群冰晶粒子的凝華增長模擬。根據Chen and Lamb（1994a）提出的冰晶粒子凝華增長理論，我們成功建立單個冰晶粒子的凝華增長模式。對于單個冰晶粒子不同溫度下的凝華增長，模擬效果較好，能夠抓住實驗觀測下所總結的冰晶粒子在不同溫度的增長變化規律（Libbrecht，2005；盛裴軒等，2003）。但與風洞實驗結果仍有一定差異，我們推測引起這些差異的主要原因有以下三個原因：（1）風洞實驗的初始粒子與模式所給的初始粒子存在差異。因為風洞實驗中初始粒子的大小不可控，其初始粒子為一系列不同的初始值，而初始粒子半徑對于冰晶增長起著重要的作用（Sheridan et al., 2009），因此這會導致模擬結果與觀測結果的差異；（2）內在增長比數據的不確定性。雖然內在增長比是一個只與溫度有關的參數，但實驗發現，內在增長比仍然受到通風因素的影響（Lamb and Scott, 1792），并且在本研究中，所采用的內在增長比來自Chen and Lamb（1994a）擬合多次不同實驗的觀測結果而得到，因此可能會導致內在增長比的不同；（3）模式中所采用的凝華密度比風洞實驗中的凝華密度大。凝華密度的偏高會導致模擬冰晶粒子質量結果偏大（Sulia and Harrington, 2011）。雖然單個冰晶粒子凝華增長模式基于Chen and Lamb（1994a）所提出的理論框架，模式結果與Chen and Lamb（1994a）單個冰晶粒子模擬結果相比，片狀冰晶的模擬結果與Takahashi and Fukuta（1988）和 Takahashi et al.（1991）觀測數據更接近，但柱狀冰晶模擬整體偏大，但由于觀測資料較少，因此很難比較模擬結果。

在群粒子的凝華增長試驗中，成功將冰晶粒子的形狀因子加入僅包括凝華增長過程的云分檔試驗模型中。在試驗模型中通過增加冰晶粒子縱橫比維度，將冰晶粒子數濃度的維度從二維增加為三維，并且為使冰晶粒子濃度在質量與縱橫比方向同時進行平流，釆用二維歐拉正定平流算法。該計算方法可使冰晶粒子形狀成功耦合到歐拉動力框架下的云模式中，這為研究冰晶粒子形狀對云微物理過程的影響及其動力過程對冰晶凝華增長的影響提供了新的途徑。

雖然采用與Chen and Lamb（1994b）不同的數值方法，但通過三種不同初始溫度下的群粒子的凝華增長，證明該試驗模式可以成功的模擬出冰晶粒子在不同溫度下質量和縱橫比的演變，并且不同的初始溫度環境對于冰晶粒子縱橫比的影響對比十分明顯。而通過不同溫度下兩種縱橫比分檔方案的對比，可看出分檔方案對于模擬效果具有較大影響，分檔越精細，模擬效果會更好，這與Leroy et al.（2007）所做的關于氣溶膠分檔方案不同的模擬結果有些相似。分檔方案的不同，會導致的模擬結果中粒子分布存在一定差異，造成這個結果的主要原因是平流計算過程中的數值擴散，即分辨率較低的方案中模擬結果更易出現較廣的分布范圍，并且模擬結果中也出現兩種不同方案的高值區數濃度高低不同，這是因為在相同的縱橫比范圍內，高分辨率的檔數較多，數濃度會相對降低。因此，縱橫比分檔分辨率高低對于今后將冰晶形狀加入分檔云模式的模擬有一定影響，因此在選擇縱橫比分檔時，根據實際計算能力，盡可能地選擇較大的縱橫比分檔。

由于在理想試驗條件下并未出現環境溫度的大幅度變化，所以不會出現冰晶從片狀向柱狀或柱狀向片狀轉化的情況，但在實際的大氣中，由于平流傳輸的作用，會形成柱帽狀冰晶，因此，在模式中暫時無法考慮該情況。

5 結論

基于CL94理論模型可以成功的模擬出單個冰晶粒子在凝華增長過程中，冰晶粒子的演變過程。釆用歐拉二維正定平流輸送法（MPDATA），在僅考慮凝華增長的歐拉數值試驗中，模型能夠模擬出群粒子在不同溫度下的冰晶質量，形狀等的增長變化特征，這為下一步將冰晶形狀因子加入到歐拉動力框架下的云模式奠定了基礎。從而更好地模擬冰晶粒子對降水、輻射等過程的影響。

致謝 感謝Jen-Ping Chen在模擬冰晶粒子形狀過程中提供的幫助。感謝Piotr K. Smolarkiewicz在MPDATA算法中提供的幫助。

（References）

Avramov A, Harrington J Y. 2010. Influence of parameterized ice habit on simulated mixed phase Arctic clouds [J]. J. Geophys. Res., 115 (D3): D03205, doi:10.1029/2009JD012108.

Bailey M, Hallett J. 2002. Nucleation effects on the habit of vapour grown ice crystals from －18 to －42°C [J]. Quart. J. Roy. Meteor. Soc., 128 (583) : 1461–1483, doi:10.1002/qj.200212858304 .

Bailey M, Hallett J. 2004. Growth rates and habits of ice crystals between －20 and －70°C [J]. J. Atmos. Sci., 61 (5): 514–544, doi:10.1175/ 1520-0469(2004)061<0514:GRAHOI>2.0.CO;2.

Castellano N E, Avila E E, Bürgesser R E, et al. 2014. The growth of ice particles in a mixed phase environment based on laboratory observations [J]. Atmos. Res., 150: 12–20, doi:10.1016/j.atmosres.2014.07.010.

Chen J P. 1992. Numerical simulations on the redistribution of atmospheric trace chemicals through cloud processes [D]. Ph. D. dissertation, The Pennsylvania State University, 342pp.

Chen J P, Lamb D. 1994a. The theoretical basis for the parameterization of ice crystal habits: Growth by vapor deposition [J]. J. Atmos. Sci., 51 (9): 1206–1222, doi:10.1175/1520-0469(1994)051<1206:TTBFTP>2.0.CO;2.

Chen J P, Lamb D. 1994b. Simulation of cloud microphysical and chemical processes using a multicomponent framework. Part I: Description of the microphysical model [J]. J. Atmos. Sci., 51 (18): 2613–2630, doi:10. 1175/1520-0469(1994)051<2613:SOCMAC>2.0.CO;2 .

Fridlind A M, Ackerman A S, McFarquhar G, et al. 2007. Ice properties of single-layer stratocumulus during the mixed-phase Arctic cloud experiment: 2. Model results [J]. J. Geophys. Res., 112 (D24): D24202, doi:10.1029/2007JD008646.

Fukuta N. 1969. Experimental studies on the growth of small ice crystals [J]. J. Atmos. Sci. 26 (3): 522–531, doi:10.1175/1520-0469(1969)026<0522: ESOTGO>2.0.CO;2.

Fukuta N. 1980. Development of fast falling ice crystals in clouds at －10°C and its consequences in ice phase processes [C]// Proc. 8th Int. Conf. on cloud Physics. Clernont-Ferrand, France, ICCP/IAMAP, 97–10.

Fukuta N, Takahashi T. 1999. The growth of atmospheric ice crystals: A summary of findings in vertical supercooled cloud tunnel studies [J]. J. Atmos. Sci., 56 (12): 1963–1979, doi:10.1175/1520-0469(1999)056< 1963:TGOAIC>2.0.CO;2.

龔乃虎. 2007. 關于冰晶雪花融化問題的研究——進展與展望 [J]. 高原氣象, 18(3): 368–376. Gong Naihu. 2007. The studies on the melting problem of ice crystals and snowflakes—Progress and prospect [J]. Plateau Meteorology (in Chinese), 18 (3): 368–376, doi:10.3321/j.issn: 1000-0534.1999.03.013.

Harrington J Y, Meyers M P, Walko R L, et al. 1995. Parameterization of ice crystal conversion processes due to vapor deposition for mesoscale models using double-moment basis functions. Part 1: Basic formulation and parcel model results [J]. J. Atmos. Sci., 52 (23): 4344–4366, doi:10.1175/1520-0469(1995)052<4344:POICCP>2.0.CO;2.

Harrington J Y, Reisin T, Cotton W R, et al. 1999. Cloud resolving simulations of Arctic stratus. Part II: Transition-season clouds [J]. Atmos. Res., 51 (1): 45–75, doi:10.1016/S0169-8095(98)00098-2.

Harrington J Y, Sulia K, Morrison H. 2013a. A method for adaptive habit prediction in bulk microphysical models. Part I: Theoretical development [J]. J. Atmos. Sci., 70 (2): 349–364, doi:10.1175/JAS-D-12-040.1.

Harrington J Y, Sulia K, Morrison H. 2013b. A method for adaptive habit prediction in bulk microphysical models. Part II: Parcel model corroboration [J]. J. Atmos. Sci., 70 (2) : 365–376, doi:10.1175/JAS- D-12-0152.1.

Hashino T, Tripoli G J. 2007. The spectral ice habit prediction system (SHIPS). Part I: Model description and simulation of the vapor deposition process [J]. J. Atmos. Sci., 64 (7): 2210–2237, doi:10.1175/ JAS3963.1.

Hashino T, Tripoli G J. 2008. The spectral ice habit prediction system (SHIPS). Part II: Simulation of nucleation and depositional growth of polycrystals [J]. J. Atmos. Sci., 65 (10): 3071–3094, doi:10.1175/ 2008JAS2615.1.

Hashino T, Tripoli G J. 2011. The spectral ice habit prediction system (SHIPS). Part III: Description of the ice particle model and the habit-dependent aggregation model [J]. J. Atmos. Sci., 68 (6): 1125–1141, doi:10.1175/2011JAS3666.1.

黃庚, 蘇正軍, 關立友, 等. 2007. 冰雪晶碰并勾連增長的實驗與觀測分析 [J]. 應用氣象學報, 18 (4): 561–567. Huang Geng, Su Zhengjun, Guan Liyou, et al. 2007. Observation and analysis of the aggregation growth among ice–snow crystals [J]. J. Appl. Meteor. Sci., 18 (4): 561–567, doi:10.3969/j.issn.1001-7313.2007.04.018.

Lamb D, Scott W D. 1972. Linear growth rates of ice crystals grown from the vapor phase [J]. J. Cryst. Growth, 12 (1): 21–31, doi:10.1016/0022- 0248(72)90333-8.

Leroy D, Wobrock W, Flossmann A I. 2007. On the influence of the treatment of aerosol particles in different bin microphysical models: A comparison between two different schemes [J]. Atmos. Res., 85(3–4): 269–287, doi:10.1016/j.atmosres.2007.01.003.

李娟, 毛節泰. 2006. 冰晶性質對卷云輻射特征影響的模擬研究 [J]. 氣象, 32 (2): 9–13. Li Juan, Mao Jietai. 2006. Simulation about influences of ice property changes on cirrus radiative properties [J]. Meteor. Mon., 32 (2): 9–13, doi:10.3969/j.issn.1000-0526.2006.02.002.

Libbrecht K G. 2005. The physics of snow crystals [J]. Reports on Progress in Physics, 68 (4): 855–895, doi:10.1088/0034-4885/68/4/R03.

Lin Y L, Farley R D, Orville H D. 1983. Bulk parameterization of the snow field in a cloud model [J]. J. Climate Appl. Meteor., 22 (6): 1065–1092, doi:10.1175/1520-0450(1983)022<1065:BPOTSF>2.0.CO;2 .

Locatelli J D, Hobbs P V. 1974. Fall speeds and masses of solid precipitation particles [J]. J. Geophys. Res., 79 (15): 2185–2197, doi:10.1029/ JC079i015p02185.

Meyers M P, Walko R L, Harrington J Y, et al. 1997. New RAMS cloud microphysics parameterization. Part II: The two-moment scheme [J]. Atmos. Res., 45 (1): 3–39, doi:10.1016/S0169-8095(97)00018-5.

Miller T L, Young K C. 1979. A numerical simulation of ice crystal growth from the vapor phase [J]. J. Atmos. Sci., 36 (3): 458–469, doi:10.1175/1520-0469(1979)036<0458:ANSOIC>2.0.CO;2.

Mitchell D L, Zhang R Y, Pitter R L. 1990. Mass-dimensional relationships for ice particles and the influence of riming on snowfall rates [J]. J. Appl. Meteor., 29 (2): 153–164, doi:10.1175/1520-0450(1990)029<0153: MDRFIP>2.0.CO;2.

Morrison H, Grabowski W W. 2008. A novel approach for representing ice microphysics in models: Description and tests using a kinematic framework [J]. J. Atmos. Sci., 65 (5): 1528–1548, doi:10.1175/ 2007JAS2491.1.

Morrison H, Grabowski W W. 2010. An improved representation of rimed snow and conversion to graupel in a multicomponent bin microphysics scheme [J]. J. Atmos. Sci., 67 (5): 1337–1360, doi:10.1175/ 2010JAS3250.1.

Morrison H, Shupe M D, Pinto J O, et al. 2005. Possible roles of ice nucleation mode and ice nuclei depletion in the extended lifetime of Arctic mixed-phase clouds [J]. Geophys. Res. Lett., 321(18): L18801, doi:10.1029/2005GL023614.

Nelson J T. 1994. A theoretical study of ice crystal growth in the atmosphere [D]. Ph. D. dissertation, University of Washington, 183pp.

Pruppacher H R, Klett J D. 1997. Microphysics of Clouds and Precipitation [M]. Dordrecht Boston: Kluwer Academic Publishers.

Reisner J, Rasmussen R M, Bruintjes R T. 1998. Explicit forecasting of supercooled liquid in winter storms using the MM5 mesoscale model [J]. Quart. J. Roy. Meteor. Soc., 124 (548): 1071–1107, doi:10.1002/qj. 49712454804 .

Ryan B F, Wishart E R, Holroyd E W. 1974. The densities and growth rates of ice crystals between ?5°C and ?9°C [J]. J. Atmos. Sci., 31 (8): 2136– 2141, doi:10.1175/1520-0469(1974)031<2136:TDAGRO>2.0.CO; 2.

Ryan B F, Wishart E R, Shaw D. E. 1976. The growth rates and densities of ice crystals between ?3°C and ?21°C [J]. J. Atmos. Sci., 33(5): 842–850, doi:10.1175/1520-0469(1976)033<0842:TGRADO>2.0.CO;2.

Sei T, Gonda T. 1989. The growth mechanism and the habit change of ice crystals growing from the vapor phase [J]. Journal of Crystal Growth, 94 (3): 697–707, doi:10.1016/0022-0248(89)90094-8.

盛裴軒, 毛節泰, 李建國, 等. 2003. 大氣物理學 [M]. 北京: 北京大學出版社. Sheng Peixuan, Mao Jietai, Li Jianguo, et al. 2003. Atmospheric Physics (in Chinese) [M]. Beijing: Beijing University Press.

Sheridan L M, Harrington J Y, Lamb D, et al. 2009. Influence of ice crystal aspect ratio on the evolution of ice size spectra during vapor depositional growth [J]. J. Atmos. Sci., 66 (12): 3732–3743, doi:10.1175/ 2009JAS3113.1.

Smolarkiewicz P K, Grabowski W W. 1990. The multidimensional positive definite advection transport algorithm: Nonoscillatory option [J]. J. Comput. Phys., 86 (2): 355–375, doi:10.1016/0021-9991(90)90105-A.

Smolarkiewicz P K. 2006. Multidimensional positive definite advection transport algorithm: An overview [J]. Int. J. Numer. Meth. Fluids, 50 (10): 1123–1144, doi:10.1002/fld.1071.

Stephens G L, Tsay S C, Stackhouse P W, et al. 1990. The relevance of the microphysical and radiative properties of cirrus clouds to climate and climatic feedback [J]. J. Atmos. Sci., 47 (14): 1742–1754, doi:10.1175/ 1520-0469(1990)047<1742:TROTMA>2.0.CO;2.

Sulia K J, Harrington J Y. 2011. Ice aspect ratio influences on mixed-phase clouds: Impacts on phase partitioning in parcel models [J]. J. Geophys. Res., 116: D21309, doi:10.1029/2011JD016298.

Sun J M, Ariya P A, Leighton H G, et al. 2012. Modeling study of ice formation in warm-based precipitating shallow cumulus clouds [J]. J. Atmos. Sci., 69 (11): 3315–3335, doi:10.1175/JAS-D-11-0344.1.

Takahashi T, Fukuta N. 1988. Supercooled cloud tunnel studies on the growth of snow crystals between －4 °C and －20°C [J]. J. Meteor. Soc. Jpn., 66 (6): 841–855.

Takahashi T, Endoh T, Wakahama G, et al. 1991. Vapor diffusional growth of free-falling snow crystals between －3 and －23°C [J]. J. Meteor. Soc. Jpn., 69 (1): 15–30.

Thompson G, Rasmussen R M, Manning K. 2004. Explicit forecasts of winter precipitation using an improved bulk microphysics scheme. Part I: Description and sensitivity analysis [J]. Mon. Wea. Rev., 132 (2): 519– 542, doi:10.1175/1520-0493(2004)132<0519:EFOWPU>2.0.CO;2.

Thompson G, Field P R, Rasmussen R M, et al. 2008. Explicit forecasts of winter precipitation using an improved bulk microphysics scheme. Part II: Implementation of a new snow parameterization [J]. Mon. Wea. Rev., 136 (12): 5095–5115, doi:10.1175/2008MWR2387.1.

Woods C, Stoelinga M, Locatelli J. 2007. The IMPROVE-1 storm of 1–2 February 2001. Part III: Sensitivity of a mesoscale model simulation to the representation of snow particle types and testing of a bulk microphysical scheme with snow habit prediction [J]. J. Atmos. Sci., 64 (11): 3927–3948, doi:10.1175/2007JAS2239.1.

周廣強, 趙春生, 秦瑜, 等. 2006. 一個簡單的雪晶輻射參數化方案及雪晶輻射對中尺度降水的影響 [J]. 應用氣象學報, 17 (2): 129–137. Zhou Guangqiang, Zhao Chunsheng, Qin Yu, et al. 2006. A simple parameterization scheme of snow-particle radiative properties and effect on mesoscale precipitation [J]. J. Appl. Meteor. Sc. (in Chinese), 17 (2): 129–137, doi:10.3969/j.issn.1001-7313.2006.02.001.

Numerical Simulations of Ice Crystal Growth in the Water Vapor Deposition Process for Typical Ice Habits

GE Sen1, 2, SUN Jiming1, 2, and NIU Shengjie1

1,,,210044,2,,,100029

A spherical shape is normally assumed for ice particles in most cloud models. However, the shapes of ice crystals in the real atmosphere are much more complicated. Ice crystal habits will change with the variation in ambient temperature and water vapor saturation. Such habits have been proven as critical parameters impacting cloud simulations. In this study, based on the theoretical model of the deposition growth of an ice crystal, we firstly simulated the growth of a single ice crystal by water vapor deposition under temperatures from ?1°C to ?30°C. The model can capture the evolution of axis length (for prism face;for basal face), mass and aspect ratio in comparison with the data from wind tunnel observations. We further simulated the water vapor deposition growth of non-spherical ice crystals with the two- dimensional positive definite advection transport algorithm (MPDATA). Furthermore, in order to test the feasibility of applying such a treatment into the Eulerian dynamical framework, the mass growth of ice crystals under different bin resolutions for the aspect ratio was simulated. The results showed that the model using the MPDATA method can capture the evolution of ice crystals for both their mass and their aspect ratio. Even though some Lagrangian models with bin microphysics involve microphysical processes for non-spherical ice crystals, their schemes with the hybrid Lagrangian– Eulerian advection transport method cannot be applied into cloud models under the Eulerian dynamical framework, which can simulate more complicated microphysical processes and dynamical processes involved in ice particles.

Ice crystal, Vapor deposition, Aspect ratio, Multidimensional positive definite advection transport algorithm (MPDATA)

10.3878/j.issn.1006-9895.1508.15151.

1006-9895(2016)03-0617-13

P40

A

10.3878/j.issn.1006-9895.1508.15151

2015-03-10；網絡預出版日期 2015-08-07

葛森，女，1990年出生，碩士研究生，主要研究方向為云霧降水物理學。E-mail: gesen_12@nuist.edu.cn

孫繼明，E-mail: jimings@mail.iap.ac.cn

國家重點基礎研究發展計劃（973計劃）項目2014CB441403，國家自然科學基金項目41375138，中國科學院大氣物理研究所百人計劃項目Y16B015601，江蘇高校優勢學科建設工程項目PADA

Founded by National Basic Research Program of China (973 Program) (Grant 2014CB441403), National Natural Science Foundation of China (Grant 41375138), One Hundred Talents Project of The Chinese Academy of Sciences (Grant Y16B015601), A Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions

葛森，孫繼明，牛生杰. 2016. 典型非球形冰晶粒子的凝華增長數值模擬試驗 [J]. 大氣科學, 40 (3): 617?629. Ge Sen, Sun Jiming, Niu Shengjie. 2016. Numerical simulations of ice crystal growth in the water vapor deposition process for typical ice habits [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 40 (3): 617?629,