基于守恒重映射的全球非结构网格模式同化及预报初步应用

陈耀登,杨登宇

南京信息工程大学 气象灾害教育部重点实验室/气象灾害预报预警与评估协同创新中心,江苏 南京 210044

全球大气数值模式对于业务天气预报和气候预测研究具有重要作用(朱跃建,2020)。近年来,随着数值计算理论的突破以及计算机技术的提升,全球大气模式在基本方程组、数值计算方法、球面网格以及物理参数化方案等方面取得了长足进步(麻巨慧等,2011;李晓莉等,2019;Randall et al.,2019)。

为解决经纬度格点模式在南北极点奇异、谱模式并行效率偏低等问题,新一代全球模式多采用非结构球面网格,如GME(Majewski et al.,2002)、FV3(Finite-Volume Cubed-Sphere Dynamical Core;Lin,2004)、NICAM(Nonhydrostatic ICosahedral Atmospheric Model;Satoh et al.,2008)、MPAS(Model for Prediction Across Scales;Skamarock et al.,2012)、CAM-SE(Community Atmosphere Model with Spectral Element dynamical core;Dennis et al.,2012)及ICON(ICOsahedral Non-hydrostatic modelling framework;Zängl et al.,2015)等。由美国国家大气研究中心(NCAR)等开发的MPAS-A(Model for Prediction Across Scales-Atmosphere)作为一款全球非静力大气模式,采用了在球面重心Voronoi网格上基于C网格离散方案的有限体积动力框架,(Thuburn et al.,2009;Ringler et al.,2010),其不仅能实现二十面体网格所具有的全球准均匀分布形态,而且支持基于点密度函数的任意地区局部加密。自发布以来,该模式受到国内外学者的广泛关注并得到了多方面的应用研究(Landu et al.,2014;Pilon et al.,2016;刘维等,2017;高元勇等,2019;Zhao et al.,2019;Kramer et al.,2020)。

另一方面,数值模式预报准确性不仅取决于模式自身的架构,也在很大程度上依赖于初始场的质量。利用资料同化技术获得更加合理的初始场,成为MPAS-A模式发展过程中必须要考虑的问题。目前已有学者针对MPAS-A的同化进行研究并取得一定成果,Ha et al.(2017)将MPAS-A模式接入到DART(Data Assimilation Research Testbed;Anderson et al.,2009)系统进行资料同化,然而DART系统基于EnKF(Ensemble Kalman Filter)方法,它的实现需要消耗大量计算资源,不便于投入业务应用。Bullock et al.(2018)则在MPAS-A中引入牛顿松弛逼近方法,以再分析数据为目标场,使模式在积分过程中向其不断逼近,但该方法只能够利用目标场中已存在的观测信息,并且其参数设置需要一定的经验积累。因此,为改进MPAS-A模式初值,相关同化研究还有待更多的开展。

与此同时,格点统计插值系统GSI(Gridpoint Statistical Interpolation)作为美国国家环境预报中心(NCEP)的业务同化系统,能够实现全球或区域范围的多源资料快速同化(Wu et al.,2002;Kleist et al.,2009)。GSI的同化框架基于变分同化方法,通过外循环和内循环迭代来解决目标代价函数的极小化问题,能够同时对不同种类的观测资料进行同化分析。GSI系统在业务运行中不断取得发展,现已支持更加先进的集合-变分混合同化技术(Wang et al.,2013;Wu et al.,2017;Duda et al.,2019),并在相关同化应用中也发挥了重要作用(董海萍等,2017;Tong et al.,2018;张涛等,2019)。将GSI与MPAS-A进行结合,即可利用GSI系统所具备的优势,为MPAS-A提供更加合理的初始场。

但是与MPAS-A基于非结构网格体系所不同的是,GSI的分析变量定义在传统的结构化高斯网格上。本研究为了发挥业务同化系统GSI成熟的多源观测资料快速同化能力,以及MPAS-A非结构网格具备的模拟优势,基于守恒重映射方法实现非结构与结构化球面网格转换,构建了一套GSI-MPAS同化及预报框架,并进行相关测试以及同化预报的初步应用。

1 GSI-MPAS同化及预报框架简介

本文的GSI-MPAS同化及预报框架如图1所示,其中同化系统为GSI,预报模式为MPAS-A,并通过转换接口实现气象要素在MPAS-A与GSI之间的转换。具体流程为:

图1 GSI-MPAS同化及预报框架流程Fig.1 Flowchart of GSI-MPAS assimilation and prediction framework

1)MPAS-A模式输出的预报场,通过转换接口进行非结构网格到结构网格的转换,并处理为GSI所需格式后,作为背景场输入到GSI系统进行同化;

2)GSI系统同化生成的分析场,经转换接口进行结构网格到非结构网格的转换,并处理为MPAS-A所需格式后,计算水平风场、位温、水汽混合比、地表气压以及干空气密度等气象要素的分析增量并更新MPAS-A的初始场,重启MPAS-A预报。

1.1 基于守恒重映射的水平网格转换

MPAS-A与GSI之间的差异首先体现在水平网格的不一致(图2),这是建立转换接口时需要首先考虑的问题。本文采用守恒重映射方法进行水平网格结构转换,其特点是能够保持物理量在转换前后的积分性质不变,对于维持模式的能量和水汽平衡具有重要意义(Jones,1999)。

图2 MPAS-A非结构网格(a)与GSI结构网格(b)的示意Fig.2 Schematic diagrams of (a)MPAS-A unstructured mesh and (b)GSI structured mesh

守恒重映射方法要求各个格点唯一对应于具有一定面积的网格单元。对于模式物理量，其全球积分表示为：

(1)

因物理量在模式网格上离散化分布，设源网格中第个格点的量值为，面积为，上式可近似表示为：

(2)

设物理量的重映射为，则目标网格与源网格满足关系：

(3)

(4)

若该目标网格共覆盖源网格中个格点，则上式进一步表示为：

(5)

其中：物理量通过下式进行估计，

(6)

由于采用C网格离散方案,MPAS-A的水平动量场定义在网格边界并与其正交(Skamarock et al.,2012)。同时,MPAS-A框架中采用径向基函数(Bonaventura et al.,2011)将网格边界的水平动量转换为网格中心的经向风和纬向风分量,以满足垂直动量方程与参数化方案的需要。转换接口在进行风场转换时,同样使用的是网格中心的风速分量。但值得注意的是,将GSI分析结果转换至MPAS-A非结构网格时,必须将网格中心的风速增量再投影至网格边界,才能完成水平动量场的更新。这一处理方案也与Ha et al.(2017)研究中所采用的方法一致。

1.2 垂直层坐标转换

在垂直方向,MPAS-A采用了跟随地形的高度坐标(Klemp,2011):

=+()(,,)。

(7)

其各个垂直层的气压通过状态方程诊断获得(Skamarock et al.,2012)。而GSI采用与业务模式GFS一致的地形追随气压坐标(Sela,2009),在其原始方案中,模式的三维气压场直接通过地表气压及两组垂直廓线系数进行确定:

(8)

因此,必须进行必要的垂直坐标转换才能使得GSI系统内部各个模式层上的气象数据能够正确对应MPAS-A模式层上的气压和高度信息。

本研究采用反距离权重插值和线性外推方法,将MPAS-A垂直方向位于模式面半层的气压数据插值到整层,并输入到GSI系统中进行计算。这一方法避免了GSI系统采用其原始方案推导模式三维气压场。图3为不同垂直层方案对于水汽混合比的处理结果与再分析资料的对比,其中图3a为FNL分析场,图3b、图3c分别表示基于GSI原始方案及本研究方案将图3a所示变量处理至模式层后的垂直分布。从图中可见,相比GSI原始方案(图3b),本研究所用方案的处理结果在垂直方向与源数据保持一致(图3c)。与气压类似,各模式层的高度信息也通过转换接口向GSI系统传递。

图3 不同垂直层方案下沿115°E的水汽混合比的垂直分布(单位:kg·kg-1):(a)FNL分析场;(b)GSI原始方案;(c)本研究方案Fig.3 Vertical distributions of water vapor mixing ratio along 115°E under different vertical coordinate schemes (units:kg·kg-1):(a)FNL analysis field;(b)GSI original scheme;(c)scheme of this paper

2 网格设置及转换功能检验

2.1 网格设置

本研究采用MPAS-A全球准均匀网格,分辨率约30 km。根据谱截断波数的不同,GSI系统中预置了多套不同分辨率的高斯结构网格。本文从中选取了T382方案,其纬向及经向格点数分别为1 152和576,分辨率约35 km,与MPAS-A所用网格最为接近。在垂直方向,MPAS-A与GSI均设置为64层,其中MPAS-A的层顶高度为30 km。

2.2 网格转换功能检验

2.2.1 检验方法

在GSI/MPAS同化及预报框架进行应用前,需首先对转换接口进行检验。检验时选取了2016年6月19日06时(世界时,下同)的GFS 0.25°分析场产品,经MPAS-A预处理模块插值到非结构网格作为基准值,然后通过转换接口将其转换至GSI结构网格后再转换回原网格,通过以下误差范数(Ullrich et al.,2009)衡量转换结果同基准值的差异:

(9)

(10)

其中，表示物理量的全球积分：

(11)

和分别表示基准值及经过两次转换后的计算结果。表征全场大尺度的误差特征，表征全场小尺度的误差特征。

2.2.2 结果对比

图4 模式物理量经过两次守恒重映射转换后的误差范数的垂直分布(a,b,c为L1范数;d,e,f为L2范数):(a,d)纬向风;(b,e)位温;(c,f)水汽混合比Fig.4 Vertical distributions of standard error of model physical quantities after two conservative remapping transformations:(a,b,c.large-scale transformation;d,e,f.small-scale transformation):(a,d)zonal wind;(b,e)potential temperature;(c,f)water vapor mixing ratio

图5 模式物理量原始场(a,b;纬向风单位:m/s;位温单位:K)以及二阶与一阶精度守恒重映射转换结果同原始场的绝对误差百分率之间的差异(c,d;单位:%):(a,c)模式第25层纬向风;(b,d)模式第25层位温Fig.5 (a,b)Original field of model physical quantities (zonal wind with the unit of m/s and potential temperature with the unit of K) and (c,d)differences in absolute error percentage between the second-order and the first-order precision conservative remapping transformation results with the original field (units:%):(a,c)zonal wind at model level 25;(b,d)potential temperature at model level 25

在应用一阶或二阶精度的守恒重映射方法进行转换时,都需要利用预先生成的转换权重系数进行计算。在一次转换过程中,采用一阶精度方法由MPAS-A非结构网格转换至GSI结构网格用时为93.996 18 s,从GSI结构网格转换回MPAS-A非结构网格用时为47.220 25 s;二阶精度方法在这两个过程中的用时分别为119.549 10 s和65.299 86 s。从计算时效性方面来看,二阶精度方法由于考虑了更多信息,在进行转换时相比一阶方法消耗更多时间,但差异并不显著。二阶精度方法的耗时主要集中在转换权重系数的生成环节,在网格确定的情况下,只需预先生成转换权重系数后便可重复使用。因此,采用一阶或二阶方法都能够在实际运行时实现快速转换。

测试分析表明,采用二阶精度的守恒重映射方法使得各变量的转换误差在大尺度和小尺度特征上相比一阶精度都有减小,故二阶精度守恒重映射方法对于变量的估计更加准确,且耗时相比一阶精度方法并未显著增加。因此,为尽量减小转换误差对结果的影响,本研究对于各变量均采用二阶精度守恒重映射方法进行网格结构转换。

3 循环同化及预报试验

为进一步分析GSI-MPAS框架的同化及预报结果,本节继续基于30 km准均匀网格,进行了为期一周的连续循环同化及预报试验,并将结果与ERA5再分析资料进行对比。

3.1 试验设计

试验选取的时段为2016年6月19日00时—2016年6月26日00时,初始场采用分辨率为0.25°的NCEP FNL数据。控制试验于起始时刻冷启动并进行预报,称作“30 km_CTRL”。循环同化试验于2016年6月19日06时采用基于初始时刻的6 h预报场作为背景场,开始进行时间间隔为6 h的连续循环同化及预报,称作“30 km_DA”。对于00时或12时同化生成的分析场,通过MPAS-A进行24 h预报,其他时刻同化完成后只进行6 h预报,以获得下一分析时刻的背景场。

MPAS-A模式为v7.0版本,采用Grell-Freitas积云对流方案(Grell and Freitas,2014),WSM6微物理过程方案(Hong and Lim,2006)、YSU边界层方案(Hong et al.,2006)、RRTMG长波及短波方案(Iacono et al.,2008)以及Noah陆面过程方案(Chen and Dudhia,2001)等。另外,在试验期间启用了海温更新,模式每逢00时从与初始场规格一致的FNL资料中获取相应时刻的海温及海冰数据。

GSI系统采用全球三维变分分析方案,其控制变量为流函数、非平衡速度势、非平衡虚温、非平衡地表气压、标准化相对湿度以及臭氧和云凝结混合比(Kleist et al.,2009),采用官方提供的全球背景误差协方差。观测数据为GDAS全球常规观测资料,其2016年6月19日12时的观测种类及位置分布如图6所示。

图6 循环同化试验于2016年6月19日12时采用的观测资料种类(括号内数字为观测数量)及其分布Fig.6 Types of observation data (numbers in brackets refer to the number of observations) and their distribution used in the cycling assimilation experiments at 1200 UTC 19 June 2016

3.2 同化及预报结果评估

3.2.1 高空及地面分析场

各组试验在最后一个同化时刻(2016年6月26日00时)的500 hPa位势高度与温度分析场与ERA5资料的差异如图7所示。从图中可见,控制试验模拟的500 hPa温度场相比ERA5资料在中纬度地区有明显偏高,而同化试验在一定程度上减小了控制试验的温度偏差,以北半球亚欧大陆南部最为显著。图7d为该同化时刻的500 hPa温度分析增量场,可见由于北半球观测资料较多(图6),分析增量主要集中在北半球地区,使得该地区的改进相比南半球更大。从高度场上看,循环同化试验的模拟结果相比控制试验更接近ERA5资料,特别是对北半球的槽脊系统定位更加准确。图7e进一步描绘了循环同化试验及控制试验的北半球500 hPa高度场与ERA5资料的距平相关系数随模拟时间的变化。从此图可见,控制试验的高度场距平相关系数随模拟时间的延伸而不断降低,而循环同化试验由于不断吸收观测数据,其高度场距平相关系数减小幅度明显低于控制试验,并在相当长的时间内保持在较高水平。

图7 2016年6月26日00时500 hPa温度场(阴影;单位:K)和高度场(等值线;单位:gpm):(a)控制试验;(b)循环同化试验;(c)ERA5资料;(d)该时刻500 hPa温度场的分析增量;(e)7 d试验期间500 hPa高度场与ERA5资料的距平相关系数随时间的变化Fig.7 500 hPa temperature (shadings;units:K) and geopotential height (contours;units:gpm) fields at 0000 UTC 26 June 2016:(a)control experiment;(b)cycling assimilation experiment;(c)ERA5 data;(d)analysis increments of 500 hPa temperature at this time;(e)variation of anomaly correlation coefficient between 500 hPa geopotential height field and ERA5 data with time during 7 d experiment

图8为各组试验的平均海平面气压模拟结果与ERA5资料的差异。可以看出,控制试验结果中存在的偏差大值区较多(图8a),反映出其对于主要地面天气系统的模拟不够准确。而在同化试验中,偏差水平在北半球有明显降低,整体体现为微弱的正偏差,在南半球依然存在少量偏差大值区,可能仍与南半球同化的观测资料数量偏少有关。

图8 2016年6月26日00时平均海平面气压场模拟结果与ERA5资料的差异(单位:hPa):(a)控制试验;(b)循环同化试验Fig.8 Differences between simulated mean sea level pressures and ERA5 data at 0000 UTC 26 June 2016 (units:hPa):(a)control experiment;(b)cycling assimilation experiment

从整体形势场上看,GSI-MPAS同化及预报框架能够有效吸收多源观测信息,并对温压场起到合理的调整作用。

3.2.2 同化及预报结果的均方根误差

为定量分析GSI-MPAS同化及预报框架的分析场及预报场特征,将循环同化及预报试验在各个00时与12时的温度和纬向风分析场及24 h预报场分别与ERA5资料在不同区域计算均方根误差并进行时间平均(图9)。从全球整体结果上看,同化试验的均方根误差在各个气压层相比控制试验都有减小,对于控制试验误差较大的气压层,改进作用也更加明显。从不同地区结果上看,两组试验在北半球的模拟误差整体低于南半球,特别是大气中低层;而且同化试验相比控制试验在北半球带来的改进幅度也大于南半球。在热带地区,控制试验在对流层中高层预报误差较大,但由于观测资料稀少,同化试验在该地区的改进作用也最小。24 h预报场中,各变量模拟误差的整体分布特征与分析场比较一致,但预报场的误差量值相比分析场有所增加。同时,同化试验在南北半球的模拟误差差异也有进一步扩大。两组试验的均方根误差评分结果进一步表明,高空观测资料的引入对于模式三维变量场具有整体改进作用,在观测资料分布密集的北半球地区带来的改进更为显著。

图9 温度(a,d;单位:K)、纬向风(b,e;单位:m/s)和相对湿度(c,f;单位:%)的模拟结果在不同地区相比于ERA5资料的均方根误差廓线:(a,b,c)分析场;(d,e,f)24 h预报场Fig.9 Profiles of root mean square error of (a,d)temperature (units:K),(b,e)zonal wind (units:m/s) and (c,f)relative humidity (units:%) simulation results compared with ERA5 data in different regions:(a,b,c)analysis field;(d,e,f)24 h prediction field

3.2.3 我国东部梅雨期降水模拟

因循环同化试验期间正值2016年梅雨期(赵俊虎等,2018),图10比较了同化预报时段内华东地区(114°～123°E,28°～36°N)梅雨降水的纬向平均值随时间的变化,这能够反映出主要雨带的南北移动情况。其中,降水实况为自动站与CMORPH降水产品融合的逐时降水数据集。从图中可见,控制试验在初期模拟的雨带相比观测有明显北移,而同化试验在一定程度上消除了控制试验的过度空报,尽管其降水强度相比观测有所减小。随着模拟时间的推移,控制试验对于降水的预测能力出现明显降低,其模拟的降水中心与实况偏离甚远。同化试验仍然能够捕捉到6月21日及23日出现的强降水事件,虽然在降水位置及强度上与实况仍有差距,但相比控制试验改进明显,显示出同化观测资料对于改善降水预报的积极影响。

图10 中国东部梅雨带(114°～123°E,28°～36°N)的纬向平均降水随时间的变化(单位:mm/h):(a)实况;(b)控制试验;(c)循环同化试验Fig.10 Variation of zonal mean precipitation with time in Meiyu belt (28°—36°N,114°—123°E) in eastern China (units:mm/h):(a)observations;(b)control experiment;(c)cycling assimilation experiment

4 结论与讨论

为改进基于非结构计算网格的全球非静力大气模式MPAS-A的模拟性能,本研究通过基于守恒重映射方法的球面网格转换和垂直坐标转换,引入GSI全球三维变分同化系统,形成了GSI-MPAS同化及预报框架。随后,开展了一系列数值试验及循环同化预报试验,得到如下结论:

1)守恒重映射方法能够实现气象要素在MPAS-A的非结构网格与GSI的高斯结构网格之间的转换。二阶精度的守恒重映射方法相比一阶精度方法的转换误差更小,能对物理量场进行更加准确的转换。

2)经过GSI同化后,MPAS-A预报得到的动力场和热力场结果在整体上获得改进,表明多源观测资料在GSI-MPAS同化及预报框架中得到了有效应用。

3)对试验结果的初步分析表明,在观测资料分布较为密集的北半球地区,同化所得分析场带来的改进作用更加明显,其各个变量的预报结果也更为合理。同时,华东地区梅雨期的降雨预报对比也表明同化后的降水预报结果在时间和空间上均与观测更加接近。

从目前初步应用结果可见,本研究构建的GSI-MPAS同化及预报框架能够利用多源观测资料得到更加合理的初始场。但通过比较不同地区的同化及预报结果也表明,因本文同化的常规观测资料大多

数位于北半球,导致南半球及热带地区的同化改进相对较小。因此,下一步将考虑同化更多观测,如卫星资料(陈耀登等,2021)等,以减小南北半球分析结果的差异,使全球整体分析结果更加合理。另外,GSI系统的背景误差及相关资料的观测误差参数设置也需进行进一步优化。

本论文的数值计算得到了南京信息工程大学高性能计算中心的计算支持和帮助。