典型非球形冰晶粒子的凝华增长数值模拟试验

葛 森 孙继明 牛生杰



典型非球形冰晶粒子的凝华增长数值模拟试验

葛 森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后，轴数值均比初始半径大的模拟结果大，轴数值均比其他初始大的模拟结果小，而在柱状冰晶增长区，则轴的模拟结果也达到最大，轴的模拟结果反而最小；在－9°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,