1 引　言

随着中国经济的发展，人工影响天气工作的重要性日趋显现。其中为抗旱蓄水所开展的人工增雨和为保障重大活动所开展的消云减雨受到各级政府的重视和广大群众的欢迎，人工影响天气的业务工作显得尤为重要。

数值天气预报模式中微物理参数化方案的不断完善使数值模式成为人工影响天气理论研究和指导人工影响作业的重要工具。在模式模拟和人工催化方面前期已有学者开展了大量的研究工作。胡志晋等（1983）为了研究层状冷云自然降水和人工催化降水的机制建立了一维非定常数值模式，考察了宏微观物理参数和冰的凝华、凇附和繁生等过程对云中要素分布和降水的影响；毛玉华等（1993）通过在二维深对流云模式中播入一定浓度的人工冰晶胚胎研究了催化对降水、降雹的影响，结果表明对于不同强度的对流云，人工引晶仅对中等积云（云顶高为6—8 km）的增雨效果较好；洪延超（1998，1999）发展了三维双参数碘化银播撒模式，并研究了催化防雹的微物理机制；Farley等（1994）运用三维云模式研究了碘化银和惰性气体在云中的运动及对云和降水的影响。Guo 等（2006，2007）建立了三维冰粒子分档播撒模式，并运用该模式对比研究了碘化银与液态二氧化碳催化的云动力和微物理效果。

随着计算机技术的飞速发展，模式分辨率不断提高，中尺度数值预报模式越来越多地应用于人工影响天气的理论研究。史月琴等（2008）、孙晶等（2010）和高茜等（2011）运用中国气象科学研究院中尺度数值模式 MM5-CAMS（Mesoscale Model 5-Chinese Academy of Meteorological Sciences CAM）通过在显式方案中直接改变冰晶数浓度的方法进行了催化数值试验；方春刚等（2009）、何晖等，2012，2013）、刘香娥等（2016）则运用加入了碘化银与云相互作用过程的中尺度模式WRF（Weather Research and Forecasting Model）、MM5 （Mesoscale Model）对对流云、层状云进行了催化模拟；Xue等（2013）对地形云地基催化进行了百米级高分辨率的WRF-LES（Large Eddy Simulation）模拟，模式中同样耦合了碘化银催化方案，可以显式分辨催化剂粒子在大气中的湍流扩散情况。

由此可见，中尺度数值模式催化模块可以详细描述催化的物理过程，其模拟的精度和模拟能力都已有了很大的提升。目前中国多地人工影响天气办公室已在业务上引进了中尺度数值模式并开发了相关的云物理产品用以指导业务。具有代表性的如北京市人工影响天气办公室的冷云催化潜力模式识别系统，该系统基于 WRF 的预报结果，可以给出具有不同催化等级的作业潜力区（嵇磊等，2014）；中国气象局人工影响天气中心的云降水精细分析和潜势预报系统，可以提供包括由数值模式、卫星、雷达探测等获得的模式预报、观测等方面的云系信息。此外，广西（张瑞波，2005；刘丽君等，2009）、河南（周毓荃等，2001，2002）、山东（王以琳等，2002，2007，2009）、湖北（袁正腾等，2012）、宁夏气象局人工影响天气办公室等均结合本地气象特点，建立了相关的人工影响天气模式或业务平台。

从文献和相关调研发现，目前中国在人工影响天气的业务工作中并没有数值模式催化产品的相关应用，对定时、定点、定量的作业信息预报也没有相关的模式产品。北京市人工影响天气办公室基于已有的研究和业务基础，借助北京睿图快速更新多尺度分析和预报系统-短时预报子系统（RMAPS-ST）的预报结果，尝试把耦合了碘化银冷云催化模块的中尺度WRF模式应用在人工影响天气业务中，以期为人工影响天气前期作业方案的制定和不断优化提供科技支撑。

2 中尺度催化数值模式

文中主模式采用美国大气研究中心和美国国家海洋和大气管理局等联合研发的新一代中尺度气象模式WRFv3.7，该模式为动力学研究、全物理过程的天气预报、空气质量预报以及区域气候模拟等提供了一个公用的模式框架。

2.1 催化模拟背景场

模拟催化的初始场和侧边界来自北京睿图快速更新和多尺度分析预报系统-临近子系统（RMAPS-ST）的实时业务预报结果。RMAPS-ST 区域设置包含两层嵌套，文中使用第二层嵌套的数据，该嵌套包含整个华北地区，水平分辨率 3 km，网格数 550×424，垂直方向设为 50 层，区域范围如图1所示。RMAPS-ST预报分析场同化了包括由全球通信系统接收的常规资料（探空、地面和飞机报）、北京地区自动气象站观测资料、京津冀地区 7 部多普勒雷达径向风观测和华北地区 29 部雷达组网拼图数据等。由于RMAPS-ST同化了多种本地观测资料，因此该数据在本研究中较美国国家环境预报中心提供的FNL再分析数据或其他再分析数据更准确（张亦洲等，2017；刘郁珏等，2019）。

2.2 催化模拟区域及参数设置

催化模拟背景场区域设置范围较大，对于北京地区的人工影响天气作业范围而言过于宽泛，并且考虑到实际计算所用的计算资源限制和业务应用时效的需要，文中进一步将催化模拟区域设置进行调整，中心点与背景场区域中心点一致，水平分辨率仍为3 km，但水平区域范围缩小，东西为230个格点，南北为260个格点，垂直方向50层，单层区域。调整后既保证了运行时效，又能覆盖飞机最大作业区域和北京地区全部地面作业点。

表1给出了调整后人工影响催化模式基本参数信息。此外，为了进行催化效果分析和后续降水机理研究，修改了WRF模式中的注册表文件（增加碘化银比含量、播撒控制开关、播撒时间和位置参数等）；修改了模式层中双参数云微物理方案（增加新方案需要的程序代码和调用该程序需要的程序代码）；在名称列表文件中增加新催化微物理方案需要的控制选项并修改WRF的后处理程序，使模式输出的碘化银比含量能够使用，最终实现催化方案和WRF模式的完全融合。为保证计算时效和进行固定点播撒，本研究将北京地区地面火箭作业站点地理信息输入到调整后的模拟区域，以便后续催化模拟时无需重新计算实际站点在模拟区域内的位置。

2.3 催化方案 2.3.1 催化机制

文中将冷云催化模块耦合到WRF中尺度数值模式莫里森双参数云物理方案中，加入人工冰核（AgI）的预报方程，主要考虑3种碘化银粒子的成核机制，即由于布朗运动和湍流扩散而发生的接触冻结核化以及冰面过饱和时水汽在人工冰核上的凝华核化和凝结冻结核化。文中假设碘化银粒子是单分散谱，半径为0.1 μm，忽略下落末速度，平均质量m_s=2.38×10⁻¹⁴ g；一个液滴只能捕获一个活化冰核发生接触核化，忽略冰质粒与碘化银粒子的碰并及碘化银粒子的光解反活化，且所有碘化银粒子在温度低于−20℃时全部活化。分别以S_bc、S_ic、S_br、S_ir、和S_dv表示人工冰核（X_S）的汇，下标第一个字母表示成核方式（b表示布朗运动，i表示湍流扩散，d表示水汽凝华），第二个字母表示与人工冰核作用的水物质（c、r和v分别表示云水、雨水和水汽）（黄燕等，1994）。矩形空间内均匀分布的碘化银粒子初始浓度（X_S0）表示模式中人工冰核的源。

不同成核机制下人工冰核的转化率分别为：

（1）云滴与碘化银粒子的接触冻结核化（S_bc、S_ic，单位：kg/(kg·s））

${{{S}}_{{\rm{bc}}}} = \frac{{\Delta {{{X}}_{\rm{s}}}}}{{\Delta {{t}}}} = - 4{\text π} {{{D}}_{\rm{s}}}{{{X}}_{\rm{s}}}{{{N}}_{\rm{c}}}{{{R}}_{\rm{c}}} = - 9.657 \times {10^{ - 17}}{{{X}}_{\rm{s}}}{{{N}}_{\rm{c}}}{{T}}$

(1)

${{{S}}_{\rm{ic}}} = \frac{{\Delta {{{X}}_{\rm{s}}}}}{{\Delta {{t}}}} = - {\text π} {{{R}}_{\rm{c}}}^2{{{X}}_{\rm{s}}}{{{V}}_{\rm{c}}}{{E}}_{{\rm{cs}}}{{{N}}_{\rm{c}}} = - 3.142 \times {10^{ - 16}}{{{X}}_{\rm{s}}}{{{N}}_{\rm{c}}}$

(2)

（2）雨滴与人工冰核的接触冻结核化（S_br、S_ir，单位：kg/(kg·s））

$ \begin{split} {{{S}}}_{\mathrm{b}\mathrm{r}}=&-4\mathrm{{\text π} }{{{D}}}_{\mathrm{s}}{{{X}}}_{\mathrm{s}}{\int }_{0}^{\infty }\frac{\,1\,}{\,2\,}{{{D}}}_{\mathrm{r}}{{{N}}}_{0\mathrm{r}}\mathrm{e}\mathrm{x}\mathrm{p}(-{\lambda}_{\mathrm{r}}{{{D}}}_{\mathrm{r}}){\rm{d}}{{{D}}}_{\mathrm{r}}\\=&-2.724\times {10}^{-10}{{{X}}}_{\mathrm{s}}{\rho }^{-0.5}{{{q}}}_{\mathrm{r}}^{0.5}{{T}} \end{split}$

(3)

$\begin{split} {{{S}}}_{\mathrm{i}\mathrm{r}}=&-{{{X}}}_{\mathrm{s}}{{{E}}}_{\mathrm{r}\mathrm{s}}{\int }_{0}^{\infty }\frac{\,{\text π}\, }{\,4\,}{{{D}}}_{\mathrm{r}}^{2}{{{V}}}_{\mathrm{r}}\left({{{ D}_{\rm r}}}\right){{{N}}}_{\mathrm{r}}\left({{{ D}_{\rm r}}}\right){{\rm{d}}}{{{D}}}_{\mathrm{r}}\\=&-2.54{{{E}}}_{\mathrm{r}\mathrm{s}}{\rho }^{-0.375}{{{{X}}}_{\mathrm{s}}{{q}}}_{\mathrm{r}}^{0.875}\end{split} $

(4)

式中，N_c是云滴浓度，取值3×10⁸ m⁻³，T为温度，R_c和V_c分别是云滴的半径和下落末速度（R_c=10 µm，V_c=1 cm/s），E_cs和E_rs分别是云、雨滴与碘化银粒子的碰并效率（E_cs=10⁻⁴，E_rs=0.5×10⁻⁴），D_s是碘化银粒子的扩散系数（D_s=KTB，其中B= $\dfrac{[ 1+(a'd/R_{\rm s})]}{6{\text π} \mu R_{\rm s}}$ ，a'=0.9，d=0.1 µm，R_s=0.1 µm， ${ \mu}=1.81\times$ 10⁻⁵ kg/（m·s），K=1.38×10⁻²³ J/K），D_r为雨滴粒径，N_0r为雨滴谱的截距，λ为雨滴谱斜率，V_r为雨滴下落末速度。ρ为空气密度，单位kg/m³，q_r为雨滴质量浓度，单位kg/kg。

（3）水汽在人工冰核上的凝华核化（S_dv，单位：kg/（kg·s））

当 5℃ ${\text{≤}} $ $\Delta T {\text{≤}}$ 20℃，

${{{S}}_{{\rm{dv}}}} = {{{m}}_{\rm{s}}}\frac{{{\rm{d}}{{{N}}_{{\rm{aD}}}}{{(\Delta {{T}})}_{}}}}{{{{\rm d}{ t}}}} = \omega \frac{{\partial [{{{X}}_{\rm{s}}}{{{N}}_{\rm{a}}}(\Delta {{T}})]}}{{\partial {{Z}}}}\dfrac{1}{{{{N}}_{\rm{a}}}(20\;{^\circ {\rm{C}}})}$

(5)

上式中 $ {\omega \dfrac{\partial \left[{{{X}}}_{\mathrm{s}}{{{N}}}_{\mathrm{a}}(\Delta {{T}})\right]}{\partial {{Z}}}}$ ＜0时，取 ${{{S}}_{{\rm{dv}}}}$ =0；

当 $\Delta {{T}}$ ＞20℃ 时，

$ {{S}}_{\rm{dv}} = {{{m}}_{\rm{s}}}{{{N}}_{{\rm{aD}}}}(\Delta {{T}}) = {{{X}}_{\rm{s}}}\frac{{{{{N}}_{\rm{a}}}(\Delta {{T}})}}{{{{{N}}_{\rm{a}}}(20\;{^\circ {\rm{C}}})}} $

(6)

$ {{N_{\rm a}}}(\Delta {{T}})=\left\{ {\begin{array}{*{20}{l}} \!\!\!\!0&{\Delta T {\text{＜}} {5^\circ }{\rm{C}}}\\ \!\!\!\!{\begin{array}{c} \!\!\!\!\!{{10}^3}\exp \Big[ - 0.022{{(\Delta {{T}})}^2} + \\0.88(\Delta {{T}}) - 3.8 \Big]\end{array}} &{{5^\circ }{\rm{C}} {\text{≤}} \Delta T {\text{＜}} {{20}^\circ }{\rm{C}}}\\ \!\!\!\!{1.6 \times {{10}^5}}&{\Delta T {\text{≥}} {{20}^\circ }{\rm{C}}} \end{array}} \right. $

(7)

式中， ${N_{{\rm{aD}}}}(\Delta T) = \dfrac{{{X_{\rm{s}}}}}{{{m_{\rm{s}}}}}\dfrac{{{N_{\rm{a}}}(\Delta T)}}{{{N_{\rm{a}}}(20\;{^\circ {\rm C}})}}$ ，是 $\Delta T$ 温度下人工冰核的活化数，单位m⁻³，Z为高度，单位为m，ω为垂直速度，单位m/s，ΔT为温度，ΔT=273.15−T。

碘化银粒子的成核作用使水成物和冰晶浓度发生变化，相应的是云滴冻结成冰晶，雨滴冻结成霰，水汽在人工冰核上凝华或凝结冻结成冰晶过程。理论上人工冰核碘化银核化除了上述3种机制还包括浸润冻结核化，但已有研究（Deshler，et al，1990；Chai，et al，1993；Li，et al，1997）显示浸润冻结核化对播撒效果影响并不显著。因此，文中所用催化模块没有考虑这种成核机制（刘香娥等，2016）。

2.3.2 用弹量计算

将计算的一定体积内过冷云水全部有效转化成较大冰相粒子（雪晶、霰粒子，粒径300 µm，密度为500 kg/m³）所需要的火箭弹数量作为模式计算的用弹量结果。计算中有如下的假设：根据WRF及MM5等中尺度模式的云方案中的设置，将云滴数浓度定为常数2.5×10⁸m⁻³；依据史月琴等（2008），何晖等（2013）计算结果，碘化银催化起效时间持续约120 min；根据美国标准大气密度与高度的关系表（Seinfeld，et al，2006）取火箭能达到的5 km高度处的标准密度为0.73 kg/m³；催化剂扩散时假设风对催化剂与云系的影响相同，即催化剂相对于云系的移动速度为0，因此只考虑催化剂的自身扩散，不考虑由于风的作用引起的催化剂平移。由于地面作业站点较多的延庆海拔在500 m左右，实际火箭作业角度为55°—65°，选择500 m海拔弹道表中的65°发射角所对应的参数作为计算用弹量的基本参数进行输入，用弹量（M）的计算公式为

${M} = \frac{{\rho \times {V} \times {\rm{LWC}}}} {{{{{m}}_q} \times {{n}} \times {{f}}(\Delta {{T}})}}$

(8)

式中， $\rho $ 为空气密度（此处取0.73 kg/m³），V为碘化银的扩散体积，根据相关文献计算火箭播撒后2 h碘化银的扩散体积为2.806 km³（周毓荃等，2014；申亿铭，1994），LWC为模式计算的云系过冷云含水量，单位为g/kg，m_q为最小雨滴质量，单位为g；n为实际作业时每枚火箭的碘化银含量，取10 g/枚； $f(\Delta T)$ 为碘化银的成核率，单位：个/g，成核率是温度的函数：

$f(\Delta T) = \frac{{{N_{\rm a}}(\Delta T)}}{{{N_{\rm a}}(20\;{^\circ {\rm C}})}} \times \frac{1}{{2.38 \times {{10}^{ - 14}}}}$

(9)

式中，ΔT=273.15−T（黄燕等，1994）。

将上述参数代入式（8）即可得到由模式计算的云水含量所需要的火箭弹数量。

2.3.3 催化方式

由2.2节可知北京地区地面火箭作业站点的经纬度信息已转换为人工影响天气催化模拟区域内的固定格点位置，在催化模拟时根据已知站点的作业信息，可直接在模式相对应格点处相应的高度、对应的作业时间、作业时长和作业剂量直接进行催化，实现固定点的静态播撒催化试验。

3 催化模拟运行

图2是催化模式基本流程示意，主要包括以下步骤：（1）运行对照试验：对照试验的初始场来自经过数据格式处理可为WRF-WPS读入的RMAPS-ST预报结果数据，按照2.3节的催化方案设置运行对照试验。（2）运行预报催化模拟：根据前述假设和参数设置，对对照试验的模拟结果进行处理，分析所有地面作业站点上空云系的微物理量，选出云水含量丰富、云系发展旺盛，具有催化潜力并持续一定时间的站点进行最优预报催化模拟。（3）实际催化模拟：将实际作业站点的信息，如作业时间、作业剂量、持续时间等输入催化模式，重新运行模拟。最后，对实际的催化模拟结果与对照试验结果进行比对，评估作业效果。同时该结果也可以评估前期的预报催化方案，给出改进和优化的意见，以提高预报催化的算法和实际作业的有效性。

催化模式运行预报催化模拟时按如下原则对站点进行筛选：（1）在站点上空火箭能达到的高度（4756—5519 m），温度低于碘化银核化临界温度（−5℃）所对应各层格点，根据该格点处的瞬时风向（若为单一风向，则取该格点及其下风向的邻近格点共2个点，若非单一风向，则取该格点及U、V下风向邻近各一个格点及合成风方向上邻近一个格点共4个点）计算云水比质量与雨水比质量之和，根据人工影响天气实际和理论作业经验，假设云中平均单格点过冷云水含量不小于0.1 g/kg则视为适合冷云催化作业；（2）当站点上空云水总含量满足条件（1），并且持续时间不少于30 min，则此站点视为优先选择；当多个站点同时满足条件（2）时，按过冷云水含量和持续时间从大到小进行优先度排序。当同一站点在多个时间段内满足条件（2）时，取间隔超过1 h的时段进行催化模拟。另外，对条件（1）中的云水、雨水含量之和运用3点后向平滑的方法进行处理，避免出现不合理的催化开始时间。以上假设考虑了云系的云微物理特征和适合催化的基本条件，既保证了催化的科学性、有效性，又考虑了模式的运行和预报时效。将筛选出的站点作业时间和作业剂量输入催化模式，进行预报催化模拟。对照和预催化模拟运行完毕后，分析模拟结果，对二者的降水信息进行比对，给出预催化效果的评估，如增减雨量、催化剂的扩散及分布等，也可以进行简单的催化机理分析。

4 个例应用及机制分析

北京人工影响天气办公室将上述算法和步骤以及计算结果初步实现了业务界面显示。系统界面可以展示北京地区的地面作业站点位置，也可以导入模式催化及其他输出结果，实现人机交互显示。系统平台还可以展示按作业条件优先度进行排序的所有作业站点信息（站点名称、作业时间、作业剂量、催化时长等），并给出所筛选站点相应催化信息在预催化模式输入名表（namelist.input）中的参数设置，实际运行中将该文件导出并输入模式，即可以进行预催化运行。

以2016年10月20日的人工影响天气作业过程为例，运用前述算法计算所有站点上空云水含量及其随时间的演变。图3为本次过程预催化时作业潜力最大、实际也进行了催化作业的白草洼站过冷云水含量演变，可见该站点20日午后及夜间云中过冷云水含量非常丰富，最大过冷云水含量达到3.2 g/kg。图4为预催化与实际催化所播撒的碘化银随系统移动在云中的扩散分布。可见预催化所用的碘化银含量要高于实际催化作业用量，催化剂在风场的作用下均向播撒的下游即北京东北部移动扩散。图5和6为模式预催化后的效果评估，以催化与对照试验10 min净降雨量之差及累计降雨量之差随时间的演变来展示，正值表示催化有增雨效果。图中可以看到预催化后约70 min，地面10 min净降雨量开始增加，10 min增雨量最大超过20×10⁶ kg；相应的地面累计降雨量也不断增加，最高达到230×10⁶ kg，预催化作业增雨效果明显。图7和8为将当日北京地面站点实际作业量进行输入的实际作业催化效果，可以看到在实际作业后，20日白天至夜间10 min净降雨量值低于10×10⁶ kg，地面累计净降雨量先增后减最高仅为45×10⁶ kg，增雨效果一般。对比前期预催化与实际作业催化的结果，预催化对20日过冷云系的开发和催化利用效果相对更好。

为分析增雨机制及预催化与实际催化效果的差异来源，此处选取地面有明显增雨效果的时段进行分析。图9中等值线为10月20日预催化与实际催化作业情况下19时10分（北京时，下同）北京东北部地面10 min增雨量的分布，可见此时二者的催化效果均为增雨，预催化时地面10 min增雨中心最大为0.031 mm，实际催化时地面10 min增雨中心最大为0.0021 mm，后续时次效果均类似，即预催化的地面增雨效果优于实际催化的增雨效果（图略）。为分析机制，给出沿图9中40.9°N云中水成物及碘化银比质量的纬向剖面分布（图10）。可以看到在该时次，云中4—6 km（温度−5—−15℃）内过冷云水含量比较丰富（0.05—0.25 g/kg），预催化在117.3°E附近的上空云中碘化银的含量要高于实际催化（图10a、b）。碘化银在丰富的过冷云水环境中活化为大量人工冰核，通过贝吉龙过程产生大量冰晶，进而产生雪晶，对应在图10c、d中117.3°E上空预催化雪的含量远高于实际催化，高空的雪下落融化增加了地面的降雨。表明在过冷水含量丰富的云层中应播撒足够量的碘化银方可实现对播撒区域内云水和水汽的充分开发，实现充分转换，提高增雨效果。这也是该个例中预催化比实际催化作业效果更好的主要原因。

5 总结与讨论

文中着眼北京地区人工影响天气业务对模式预报的需求，结合地基作业的特点，进行了地基碘化银冷云催化数值模式在业务应用中的相关研究，结果如下：

（1）催化模拟的背景场来自同化了多种观测资料的RMAPS-ST预报结果数据。模拟的区域范围覆盖北京地区所有的地面作业站点以及华北人工影响天气飞机作业区域。碘化银冷云催化模块耦合在WRF Morrison双参数微物理方案中，主要考虑了冰核的接触冻结核化以及水汽在人工冰核上的凝华核化。

（2）模拟催化流程主要包括：（a）运行对照试验，根据算法挑选出适合的作业站点和作业时间、作业量；（b）根据（a）给出的作业信息运行预催化模拟，给出预催化的作业效果（增、减雨量）；（c）将实际作业信息输入进行实际作业重催化模拟，同样给出作业效果。对比预催化和实际催化的效果，评估预催化方案。

需要指出的是，文中只是对冷云催化数值模式在业务中的应用进行了初步的研究。文中的理论计算如用弹量等的计算过程中有诸多假设，再加上实际的人工影响天气作业受多种主、客观因素的限制，如空域申请等。因此，北京地区的冷云催化数值模式业务应用目前在实施中有诸多不足，通过对个例的应用研究和分析预催化与实际催化累计降雨量差产生的原因可以对算法指标进行评估改进，以更好地指导应用，今后需要进一步改进和完善。此外，成熟的业务应用要涉及到诸如计算时效、预报时效、地面（或空中）作业资源配置等多种因素，如何兼顾多种影响因素并在有限的预报时段内给出合理的催化方案还需要长时间序列的个例检验。针对不同的天气形势、不同季节、不同类型的降水等催化方案都会有一定的差异，今后还需要更多的科学研究和业务应用来进行总结归纳和优化提升。