基于四维变分同化的“6·23”阜宁龙卷大涡模拟研究

刘达之 姚聃,2 梁旭东

(1.中国气象科学研究院灾害天气国家重点实验室，北京 100081; 2.中国气象局气象探测中心，北京 100081)

引言

龙卷是破坏力极强的小尺度系统，其发生往往造成巨大的人员伤亡和财产损失，产生重大社会影响[1-5]。龙卷的预报预警是气象科研和业务中急需解决的关键性难题[6]。龙卷预警主要依靠多普勒雷达对于中气旋和龙卷涡旋特征(Tornadic Vortex Signature，TVS)的识别。2016年，中国气象局启动了龙卷预警业务体系建设试验[7]，近几年该试验取得了初步成效[8]，并已在2021年5月14日江苏苏州、湖北武汉两例EF3级强龙卷等重大灾害性天气过程的预警中发挥了重要作用。广东佛山市龙卷风研究中心也开展了针对龙卷的预警业务[9-10]。然而，基于雷达监测的预警在龙卷母体风暴生成和发展以后才能发现，预警时效较短，难以满足龙卷灾害防御的需求。

龙卷预报主要依赖于数值模式对于龙卷母体风暴及其所产生环境条件的模拟和预报。数值模拟对于提高龙卷的预报时效具有重要意义，但相关研究难度较大[11]，主要原因有3点：①针对龙卷尺度的模拟技术和高时空分辨率所需运算资源的限制;②针对龙卷复杂结构和发生机理的认知不足；③龙卷过程可预报性的局限。强对流系统的触发在数值模拟和预报中是难点，尺度更小的龙卷过程和机制则较为复杂，需要以典型个例作为切入点，开展深入的系统研究。Hanley等[12]以造成严重伤亡的2013年5月20日美国俄克拉荷马州Moore市EF5级龙卷为切入点，使用100 m分辨率成功模拟出龙卷超级单体以及近地面强风，但资料同化的缺少造成位置和时间偏差明显。Snook等[13]利用集合卡尔曼滤波(Ensemble Kalman Filter，EnKF)同化雷达资料，成功通过大涡模拟再现出了50 m分辨率的龙卷结构，并基于集合预报试验分析了龙卷路径和强度对于模式的高度敏感性。Zhang等[14-15]同样以EnKF为基础，研究模式微小差异在风暴发展期间的非线性增长，并从实际和可预报性的角度分析了龙卷超级单体数值模拟和预报所面临的核心问题。以上研究为龙卷的数值模拟和预报可行性提供了有利参考。

2016年6月23日，江苏盐城发生强度达到EF4级的超强龙卷过程，在阜宁县造成98人死亡[2]，是该地区40多年以来最严重的一次龙卷灾害[16]。针对该龙卷过程，学者们采用不同方法开展了数值模拟研究。Yao等[17]以WRF模式所模拟的环境场提取垂直廓线生成理想化初始场，驱动CM1模式开展25 m分辨率的大涡模拟，研究精细化的龙卷漏斗云三维结构及其演变过程。Sun等[18]采用WRF模式进行龙卷模拟对水平分辨率的依赖性研究，指出在49 m网格上通过大涡模拟可以解析龙卷结构。该结果为进一步开展的龙卷结构和发生机理研究提供了参考，但其较大的计算资源消耗难以应用于业务预报，同时龙卷发生位置和时间存在偏差。以雷达资料同化为基础，陈锋等[19]利用雷达反射率以及径向风资料，对比同化不同资料的结果，发现同化径向速度的方案能显著改进阜宁附近涡旋的发生、发展过程模拟。李佳等[20]利用华东区域模式快速循环同化系统，对阜宁龙卷天气过程进行模拟，预报时效为30 min，模拟出具有一定强度的对流单体。但这两个研究仅使用3 km分辨率进行模拟，并没有分析在此环境场能否产生龙卷。

以上研究表明，龙卷的数值模拟既需要资料同化的改进以减小误差，也需要百米量级分辨率以解析龙卷尺度结构，龙卷的数值预报则更需要在此基础上兼顾运算效率。本研究基于WRF模式与四维变分同化方案(4DVar)系统，同化高时空分辨率的地面观测和雷达径向风资料，开展龙卷可分辨尺度大涡模拟，分析阜宁龙卷及其母体风暴的发生发展过程，并探讨基于粗分辨率判断龙卷发生条件的可能，以期为龙卷预报提供参考。

1 资料与方法

1.1 资料来源

模式初始场选用了NCEP/GFS再分析资料2016年6月23日03:00 UTC预报场，以开展2016年6月23日阜宁龙卷的模拟。同化所使用的资料包括全国地面2063站逐小时观测(气压/海平面气压、气温、露点温度、2 min平均风向和风速)以及14部多普勒天气雷达(盐城、连云港、淮安、上海青浦、南京、南通、常州、泰州、宿迁、临沂、合肥、蚌埠、马鞍山、铜陵)。天气形势分析选用NCEP/FNL最终分析资料。累积降水和地面风场分析使用了江苏省内地面分钟级观测资料。雷达观测分析使用了江苏、淮安、连云港3部雷达进行资料合成。

1.2 试验设计

采用WRF模式3.9版本，最外层网格的中心点位置为33°N、115°E，水平分辨率为9 km。模拟分为基础试验和加密试验，基础试验向内依次嵌套3 km，1 km和333 m网格，共计四重。加密试验分别针对模拟出的两个对流单体继续嵌套111 m最内层网格，共计五重(图1)。垂直方向共设置51层，其中最低层高度约为33 m。

在物理过程参数化方面，基于对前人研究[18-20]的综合测试比较，选取了Morrison双参云物理方案、RRTMG短波和长波辐射方案、Pleim-Xiu 近地面层方案和陆面过程计算方案、ACM2边界层方案以及SGS湍流参数化方案组合，关闭了积云参数化。在333 m和111 m分辨率区域采用大涡模拟。

D01～D04区域水平分辨率依次为9 km、3 km、1 km和333 m，D05、D06分别为对单体1和单体2进行加密的区域水平分辨率111 m的最内层网格图1 模式区域设置Fig.1 Mode region settings

采用WRF模式的4DVar系统开展同化试验。在03:00—04:00 UTC同化窗内，将初始场与全国地面逐小时观测资料以及周边区域雷达反演风场数据融合。雷达径向速度的处理采用Liang等[21]的方法，反演产品水平分辨率为0.1°、垂直分辨率为500 m、垂直方向为15层、时间间隔为12 min的三维风场资料。

1.3 诊断参量

龙卷的模拟中，两个常用的垂直积分参量容易混淆，表征龙卷生成环境条件的风暴相对螺旋度(Storm Relative Helicity，SRH)和表征风暴单体中旋转强度的上升气流螺旋度(Updraft Helicity，UH)。SRH为指定垂直层次内(通常取0～3 km)水平风的垂直切变及其旋转程度[22]，最初用于龙卷临近探空曲线的诊断，也可以用于判定分析模式场中有利于龙卷生成的范围。UH计算的则是指定垂直层次内(通常选取2～5 km)涡度与垂直速度的乘积[23]，用于表征三维风场中已经形成的旋转上升气流的强度，适用于模式分辨率难以解析龙卷的情况下，判定中气旋及龙卷可能发生的位置和路径。

此外，对于龙卷可分辨尺度模拟而言，类龙卷涡旋可以直接作为龙卷近地面环流的表征。参照Schenkman等[24]研究，本文采用的类龙卷涡旋标准为：低层存在闭合环流，持续时间不低于2 min，中心最大垂直涡度大于0.2 s-1，且水平风速不小于29 m·s-1。

2 结果分析

2.1 阜宁龙卷概况

2016年6月23日阜宁龙卷发生在中国东部的梅雨期，其天气背景有利于暴雨发生[25]。地面暖锋南侧高温、高湿的不稳定气层提供了良好的热力条件，高低空急流的耦合提供了良好的动力配置。由NCEP/FNL资料可知(图2a和图2b)，伴随东北冷涡的槽线南端位于渤海湾。阜宁在西太平洋副热带高压北侧，低空受西南暖湿气流影响，湿度达80%。而在对流层中上层，江苏北部受东北气流控制，高低空不同性质的气流加剧了层结不稳定。这些有利环境条件在资料同化过程中均得到了改善(图2c和图2d)。

阴影区为850 hPa湿度，单位为%；黑色实线为500 hPa位势高度，单位为gpm；红色虚线为500 hPa温度，单位为℃；箭头为850 hPa风场；蓝色三角为盐城阜宁县位置图2 2016年6月23日 00:00 UTC(a)、06:00 UTC(b)FNL最终分析资料、以及同化前(c)与同化后(d)06:00 UTC高空形势场模拟结果Fig.2 FNL final analysis data on 0000 UTC (a) and 0600 UTC (b),and the high-level field on 0600 UTC in pre-assimilation (c) and post-assimilation (d) simulation results on June 23,2016

此次阜宁龙卷发生于超级单体中，该超级单体位于西北—东南走向的线状中尺度对流系统南端[4]，与2012年7月21日北京通州张家湾龙卷相似[26]。观测资料表明，05:00—06:00 UTC时段降水达到最强，分南北两个降水极值中心，偏南侧的降水中心影响阜宁，1 h累积降水量达到25 mm以上；06:00—07:00 UTC降水区仍然影响阜宁，但强度有所减弱。产生龙卷的超级单体并未发生在最强降水区，但地面风场有明显的辐合中心。

产生龙卷的超级单体在发生发展过程中先后与多个对流单体合并，其低层有钩状回波和入流缺口，可识别出中气旋和龙卷涡旋特征(TVS)[27]。根据盐城、连云港、淮安雷达合成的1.5 km高度反射率因子，06:00—06:36 UTC，对流系统整体向东北方向移动，区域内具有南北两个强度达50 dBz的对流单体，南侧对流单体即为产生龙卷的超级单体。

2.2 龙卷超级单体模拟结果

观测资料同化前，2016年6月23日各高度背景场为均匀的偏南风。观测资料同化后，分析场中出现了明显的中尺度涡旋场。在分析场850 hPa上，山东南部、江苏北部、渤海湾出现闭合气旋性风场。江苏境内各高度偏北风分量明显增强(图略)，这使得北方冷空气进一步南压，与南方暖湿气流交汇，促进对流系统发生发展。降水方面，对比同化前后初始场预报6 h累计降水可知，同化前，在阜宁附近没有降水，而同化后出现了降水集中区域，与观测更为接近(图3)。

阴影区为6 h累积降水量；红框标记所关注区域；黑点为阜宁龙卷发生位置图3 2016年6月23日03:00—09:00 UTC同化前(a)和同化后(b)模拟的累积降水量及阜宁地区累积降水量观测实况(c)Fig.3 Pre-assimilation (a) and post-assimilation (b) simulated accumulated precipitation at 0300-0900 UTC and observations in Funing region (c) on June 23,2016

图3红色方框内，1.5 km高度上06:30 UTC雷达反射率因子表明，有两个对流单体，分别命名为单体1和单体2(图4)。06:30 UTC，单体2影响阜宁，强度为40 dBz，在时间和空间上与观测值十分接近。单体1位于阜宁以西，最大回波强度略大于单体2。加密试验中，尽管单体2是较强的对流单体，但没有出现类龙卷涡旋。单体1生成了类龙卷涡旋，并出现明显的钩状回波、风场辐合和强的上升气流螺旋度(UH)(图5)。基础试验中，单体1的雷达反射率因子与UH随时间不断增强。08:30 UTC，虽然滞后于龙卷发生的实际时间，但该单体产生了较强的UH并伴有低层风的气旋性辐合，风速大于20 m·s-1。在形态上，单体1的形态呈钩状回波特征。

阴影区为1.5 km高度雷达反射率因子；红色方框分别标记单体1和单体2；黑点为阜宁龙卷发生位置图4 2016年6月23日06：30 UTC同化后模拟的阜宁地区雷达反射率因子Fig.4 Radar reflectance factors at 0630 UTC in post-assimilation simulation in Funning region on June 23,2016

加密试验的大涡模拟表明，08:30 UTC该单体具有明显的钩状回波结构，并伴有UH大值区(图6a和图6c)，说明存在强烈旋转上升气流。UH大于400 m2·s-2的极值区沿着钩状回波尾端外侧的辐合线呈环状分布，最大值已超过1200 m2·s-2。08:28—08:34 UTC,200 m高度垂直涡度以及模式最低层风场逐2 min时间演变表明(图7)，这个时间段内，涡度场分布呈圆环状，08:30 UTC涡度最大值达0.2 s-1并维持4 min，08:30—08:34 UTC涡度极值区周边的风速达40 m·s-1，并形成闭合环流，符合类龙卷涡旋特征。

单体2在基础试验和加密试验中，单体2发展最强的时段(06:20—06:40 UTC)内均未出现钩状回波、大于50 m2·s-2的UH或风场的气旋性辐合，无类龙卷涡旋生成。

在阜宁龙卷的已有研究中，Sun等[18]采用未同化观测资料模拟出了龙卷涡旋结构，07:12 UTC该龙卷涡度中心在144 m水平分辨率上，位于34°11.5′N、118°45′E，与阜宁龙卷相比，时间偏晚、位置偏北。在李佳等[19]、陈锋等[20]采用资料同化的研究中，阜宁龙卷周边区域，在相近的时间出现了较强的对流系统，但没有模拟出类龙卷涡旋。本文中，资料同化后得到了与阜宁龙卷发生位置和时间相近的较强对流单体(单体2)，但对其采用高分辨率加密不能模拟出类龙卷涡旋。而在西侧区域较晚的时间，在对单体1加密后模拟出了类龙卷涡旋。

因此，通过资料同化改进初始场，形成有利于龙卷发生的环境场，有利于模拟出类龙卷涡旋。对于龙卷发生的具体时间和位置，仍具有较大的不确定性。对于以NCEP/GFS预报场为基础，尝试对龙卷进行模拟和预报的目标而言，通过资料同化可以提高对流触发及龙卷预报的技巧，也表明模拟结果与预报业务目标存在较大差距。

阴影区为1.5 km高度雷达反射率因子；箭头为最低层风场，等值线UH=800 m2·s-2图5 2016年6月23日08:20UTC (a)、08:30 UTC (b)、08:40 UTC (c) 阜宁地区333 m分辨率网格单体1对应的雷达反射率因子、最低层风场及上升气流螺旋度Fig.5 Radar reflectivity factor,wind at the first model level,and updraft helicity of convective cell 1 at 0820 UTC (a),0830 UTC (b) and 0840 UTC (c) on the 333 m grid on June 23,2016 of Funing region

图a和图b阴影区为1.5 km高度雷达反射率因子；图c和图d阴影区为上升气流螺旋度；箭头为最低层风场；等值线为40 dBz雷达反射率因子；图c和图d分别为图a和图b中黑框范围图6 2016年6月23日阜宁地区111 m分辨率网格双向嵌套(a)、单向嵌套(b)单体1对应的雷达反射率因子、图a黑框区域(c)、图b黑框区域(d)最低层风场及上升气流螺旋度Fig.6 Radar reflectivity factor,lowest layer wind and updraft helicity of convective cell 1 in two-way (a) ,one-way(b),black box area of figure 6a (c) and figure 6b (d) nesting on the 111 m grid on June 23,2016 of Funing region

2.3 影响龙卷生成的多尺度条件

2.3.1 有利龙卷生成的环境条件

不稳定能量方面，在模式的起报初始时间(2016年6月23日03:00 UTC)，包括盐城阜宁区域在内，江苏大部分区域CAPE值达2000 J·kg-1以上(图略)，超过了龙卷性超级单体发生条件的中位数(约1800 J·kg-1)[28]，表明环境场具备充足的不稳定能量。CAPE是龙卷性超级单体生成的基本前提条件，虽然阜宁当日大部分区域CAPE值均超过龙卷性超级单体发生的中位数，此次阜宁龙卷并未发生在CAPE值最大的区域。

水平风矢量的垂直切变方面，单体1和单体2所在区域的环境，0～6 km切变均已超过20 m·s-1，而0～1 km切变则差异明显(图略)。据Thompson等[28]统计，0～6 km垂直风切变在超级单体与非超级单体间差异显著，超级单体0～6 km切变中位数约为22 m·s-1；而强龙卷、弱龙卷、非龙卷性超级单体以及非超级单体对流的0～1 km切变中位数分别为9.8 m·s-1、8.1 m·s-1、6.4 m·s-1和2.9 m·s-1。本研究中，单体1的0～1 km垂直风切变在05:30—07:00 UTC为5～7 m·s-1，在07:10—08:00 UTC为8～15 m·s-1，而单体2在05:30—06:00 UTC仅为1～2 m·s-1。即单体1的环境有利于龙卷超级单体生成，而单体2的环境更有利于非超级单体对流的生成。两个对流单体0～6 km垂直风切变均有利于超级单体发生，但仅在单体1及其周边范围内，0～1 km垂直风切变有利于龙卷的发生。

相对于低层垂直风切变而言，风暴相对螺旋度(SRH)则更集中在单体1附近范围内，呈现出区别于周边区域有利于龙卷生成的环境(图8)。 单体1的SRH达100 m2·s-2以上，单体2的SRH不足25 m2·s-2。根据Thompson等[28]统计，强龙卷、弱龙卷和非龙卷超级单体的SRH中位数依次为223 m2·s-2、184 m2·s-2和146 m2·s-2。此次阜宁龙卷SRH随模拟分辨率提升而提高，但整体分布不变。3 km分辨率中单体1的SRH07:50 UTC超过100 m2·s-2，1 km分辨率达200 m2·s-2，333 m分辨率达225 m2·s-2，但是空间分布没有发生变化。与之相比，单体2的SRH整个阶段在各分辨率区域均未达到25 m2·s-2。

阴影区为0～3 km SRH；三角形(正方形)为单体1(单体2)达到回波强度最大值位置图8 2016年6月23日07:50 UTC 3 km(a)、1 km(b)、333 m(c)、111 m(d)、08:00 UTC 3 km(e)、1 km(f)、333 m(g)、111 m(h)、08:10 UTC 3 km(i)、1 km(j)、333 m(k)、111 m(l)分辨率网格阜宁地区风暴相对螺旋度变化Fig.8 Temporal variations in SRH at 0750UTC in grid resolution of 3 km (a),1 km (b),333 m (c) and 111 m (d),0800 UTC in grid resolution of 3 km (e),1 km (f),333 m (g) and 111 m (h),0810 UTC in grid resolution of 3 km (i),1 km (j),333 m (k) and 111 m (l) in Funing region on June 23,2016

此次阜宁龙卷模拟，单体1和单体2前期为强度相似的两个对流单体。即使在3 km分辨率下，也已表明其中一个单体所在环境场更有利于龙卷的发生。这种环境场的差异，单体1在高分辨率模拟区域中心出现了类龙卷涡旋。模式分辨率的提升，使得SRH的量值增大，但在分布上没有明显差异。表明SRH对于龙卷的生成具有指示性，且在不使用高分辨率的条件下，就可以通过环境场的模拟分析，为能否生成龙卷提供预报依据。这体现了基于数值模式开展龙卷预报(如潜势联合概率预报[29])的可能。

2.3.2 小尺度系统的反馈作用

具备了有利的环境条件，龙卷不一定发生，且其中所涉及的过程较为复杂。为说明这一现象，通过对比测试，分析在完全相同的环境条件下，单向嵌套网格与双向嵌套网格所产生模拟结果之间的差异。即在以上针对单体1加密试验的基础上，关闭内层网格对外层网格的反馈，分析高分辨率模拟产生的小尺度过程对龙卷超级单体的影响。

采用单向嵌套后，加密试验的1.5 km高度，UH最大值达600 m2·s-2(图6b和图6d)，远低于双向嵌套试验中UH最大值(图6a和图6c)。在1.5 km高度回波反射率因子形态结构上未出现明显钩状回波，未出现持续时间超2 min、强度大于0.2 s-1的垂直涡度，未达到类龙卷涡旋标准。对比1 km、333 m和111 m不同分辨率下模式最低层最大垂直涡度可知(图9)，随着分辨率的提高，垂直涡度的强度也逐步提高。总体上，双向嵌套的垂直涡度大于单向嵌套的垂直涡度。111 m分辨率下，双向嵌套下垂直涡度的最大值于08:30 UTC达到0.2 s-1，维持4 min。单向嵌套时仅在08:40 UTC后有短暂的时间段达到0.2 s-1。

图9 2016年6月23日不同分辨率双向嵌套与单向嵌套的垂直涡度最大值变化Fig.9 Temporal variations in the maximum vertical vorticity of two-way and one-way nesting on June 23,2016

以上结果表明，即使在环境条件有利于龙卷生成的情况下，类龙卷涡旋的生成也需要小尺度系统的反馈作用，甚至超级单体的维持和发展也可能受到影响。这其中可能反映了龙卷生成理论中近地面涡旋增强机制的作用[11]，以及摩擦和湍流等物理过程的影响[24]。

3 结论与讨论

(1)同化高频次的地面观测和雷达风场资料，能有效改进模式初始场，优化环境条件的模拟，从而提升强对流系统生成位置和持续时间的模拟效果，以及对生成龙卷的超级单体风暴发生、发展的模拟能力。

(2)尽管风暴相对螺旋度等环境参量值随模式分辨率有所提高，但其空间分布并无显著差异。与对流有效位能(CAPE)和垂直风切变相比，风暴相对螺旋度的高值分布更为集中，是区分相似对流单体能否产生龙卷的较好指示因子。

(3)复杂的小尺度过程对龙卷的实际发生具有重要影响。在缺少小尺度系统的反馈作用时，即使有利于龙卷生成的环境场也可能模拟不出类龙卷涡旋。这既涉及龙卷发生的复杂机理，也对龙卷预报具有重要参考意义，需深入开展研究。

(4)基于数值模式开展龙卷预报的主要难点在于，风暴单体的对流触发问题已十分复杂，在此基础上还要准确筛选出可能产生龙卷的单体。在此过程中，无论大尺度环境场、中尺度母体风暴还是小尺度龙卷的触发，均存在难以忽视的不确定性和可预报性问题[6]。开展基于集合的数值模拟研究很可能是解决龙卷模式预报难题的有效途径。