基于守恒重映射的全球非結(jié)構(gòu)網(wǎng)格模式同化及預(yù)報(bào)初步應(yīng)用

陳耀登,楊登宇

南京信息工程大學(xué) 氣象災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室/氣象災(zāi)害預(yù)報(bào)預(yù)警與評(píng)估協(xié)同創(chuàng)新中心,江蘇 南京 210044

全球大氣數(shù)值模式對(duì)于業(yè)務(wù)天氣預(yù)報(bào)和氣候預(yù)測(cè)研究具有重要作用(朱躍建,2020)。近年來,隨著數(shù)值計(jì)算理論的突破以及計(jì)算機(jī)技術(shù)的提升,全球大氣模式在基本方程組、數(shù)值計(jì)算方法、球面網(wǎng)格以及物理參數(shù)化方案等方面取得了長(zhǎng)足進(jìn)步(麻巨慧等,2011;李曉莉等,2019;Randall et al.,2019)。

為解決經(jīng)緯度格點(diǎn)模式在南北極點(diǎn)奇異、譜模式并行效率偏低等問題,新一代全球模式多采用非結(jié)構(gòu)球面網(wǎng)格,如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)等。由美國(guó)國(guó)家大氣研究中心(NCAR)等開發(fā)的MPAS-A(Model for Prediction Across Scales-Atmosphere)作為一款全球非靜力大氣模式,采用了在球面重心Voronoi網(wǎng)格上基于C網(wǎng)格離散方案的有限體積動(dòng)力框架,(Thuburn et al.,2009;Ringler et al.,2010),其不僅能實(shí)現(xiàn)二十面體網(wǎng)格所具有的全球準(zhǔn)均勻分布形態(tài),而且支持基于點(diǎn)密度函數(shù)的任意地區(qū)局部加密。自發(fā)布以來,該模式受到國(guó)內(nèi)外學(xué)者的廣泛關(guān)注并得到了多方面的應(yīng)用研究(Landu et al.,2014;Pilon et al.,2016;劉維等,2017;高元勇等,2019;Zhao et al.,2019;Kramer et al.,2020)。

另一方面,數(shù)值模式預(yù)報(bào)準(zhǔn)確性不僅取決于模式自身的架構(gòu),也在很大程度上依賴于初始場(chǎng)的質(zhì)量。利用資料同化技術(shù)獲得更加合理的初始場(chǎng),成為MPAS-A模式發(fā)展過程中必須要考慮的問題。目前已有學(xué)者針對(duì)MPAS-A的同化進(jìn)行研究并取得一定成果,Ha et al.(2017)將MPAS-A模式接入到DART(Data Assimilation Research Testbed;Anderson et al.,2009)系統(tǒng)進(jìn)行資料同化,然而DART系統(tǒng)基于EnKF(Ensemble Kalman Filter)方法,它的實(shí)現(xiàn)需要消耗大量計(jì)算資源,不便于投入業(yè)務(wù)應(yīng)用。Bullock et al.(2018)則在MPAS-A中引入牛頓松弛逼近方法,以再分析數(shù)據(jù)為目標(biāo)場(chǎng),使模式在積分過程中向其不斷逼近,但該方法只能夠利用目標(biāo)場(chǎng)中已存在的觀測(cè)信息,并且其參數(shù)設(shè)置需要一定的經(jīng)驗(yàn)積累。因此,為改進(jìn)MPAS-A模式初值,相關(guān)同化研究還有待更多的開展。

與此同時(shí),格點(diǎn)統(tǒng)計(jì)插值系統(tǒng)GSI(Gridpoint Statistical Interpolation)作為美國(guó)國(guó)家環(huán)境預(yù)報(bào)中心(NCEP)的業(yè)務(wù)同化系統(tǒng),能夠?qū)崿F(xiàn)全球或區(qū)域范圍的多源資料快速同化(Wu et al.,2002;Kleist et al.,2009)。GSI的同化框架基于變分同化方法,通過外循環(huán)和內(nèi)循環(huán)迭代來解決目標(biāo)代價(jià)函數(shù)的極小化問題,能夠同時(shí)對(duì)不同種類的觀測(cè)資料進(jìn)行同化分析。GSI系統(tǒng)在業(yè)務(wù)運(yùn)行中不斷取得發(fā)展,現(xiàn)已支持更加先進(jìn)的集合-變分混合同化技術(shù)(Wang et al.,2013;Wu et al.,2017;Duda et al.,2019),并在相關(guān)同化應(yīng)用中也發(fā)揮了重要作用(董海萍等,2017;Tong et al.,2018;張濤等,2019)。將GSI與MPAS-A進(jìn)行結(jié)合,即可利用GSI系統(tǒng)所具備的優(yōu)勢(shì),為MPAS-A提供更加合理的初始場(chǎng)。

但是與MPAS-A基于非結(jié)構(gòu)網(wǎng)格體系所不同的是,GSI的分析變量定義在傳統(tǒng)的結(jié)構(gòu)化高斯網(wǎng)格上。本研究為了發(fā)揮業(yè)務(wù)同化系統(tǒng)GSI成熟的多源觀測(cè)資料快速同化能力,以及MPAS-A非結(jié)構(gòu)網(wǎng)格具備的模擬優(yōu)勢(shì),基于守恒重映射方法實(shí)現(xiàn)非結(jié)構(gòu)與結(jié)構(gòu)化球面網(wǎng)格轉(zhuǎn)換,構(gòu)建了一套GSI-MPAS同化及預(yù)報(bào)框架,并進(jìn)行相關(guān)測(cè)試以及同化預(yù)報(bào)的初步應(yīng)用。

1 GSI-MPAS同化及預(yù)報(bào)框架簡(jiǎn)介

本文的GSI-MPAS同化及預(yù)報(bào)框架如圖1所示,其中同化系統(tǒng)為GSI,預(yù)報(bào)模式為MPAS-A,并通過轉(zhuǎn)換接口實(shí)現(xiàn)氣象要素在MPAS-A與GSI之間的轉(zhuǎn)換。具體流程為:

圖1 GSI-MPAS同化及預(yù)報(bào)框架流程Fig.1 Flowchart of GSI-MPAS assimilation and prediction framework

1)MPAS-A模式輸出的預(yù)報(bào)場(chǎng),通過轉(zhuǎn)換接口進(jìn)行非結(jié)構(gòu)網(wǎng)格到結(jié)構(gòu)網(wǎng)格的轉(zhuǎn)換,并處理為GSI所需格式后,作為背景場(chǎng)輸入到GSI系統(tǒng)進(jìn)行同化;

2)GSI系統(tǒng)同化生成的分析場(chǎng),經(jīng)轉(zhuǎn)換接口進(jìn)行結(jié)構(gòu)網(wǎng)格到非結(jié)構(gòu)網(wǎng)格的轉(zhuǎn)換,并處理為MPAS-A所需格式后,計(jì)算水平風(fēng)場(chǎng)、位溫、水汽混合比、地表氣壓以及干空氣密度等氣象要素的分析增量并更新MPAS-A的初始場(chǎng),重啟MPAS-A預(yù)報(bào)。

1.1 基于守恒重映射的水平網(wǎng)格轉(zhuǎn)換

MPAS-A與GSI之間的差異首先體現(xiàn)在水平網(wǎng)格的不一致(圖2),這是建立轉(zhuǎn)換接口時(shí)需要首先考慮的問題。本文采用守恒重映射方法進(jìn)行水平網(wǎng)格結(jié)構(gòu)轉(zhuǎn)換,其特點(diǎn)是能夠保持物理量在轉(zhuǎn)換前后的積分性質(zhì)不變,對(duì)于維持模式的能量和水汽平衡具有重要意義(Jones,1999)。

圖2 MPAS-A非結(jié)構(gòu)網(wǎng)格(a)與GSI結(jié)構(gòu)網(wǎng)格(b)的示意Fig.2 Schematic diagrams of (a)MPAS-A unstructured mesh and (b)GSI structured mesh

守恒重映射方法要求各個(gè)格點(diǎn)唯一對(duì)應(yīng)于具有一定面積的網(wǎng)格單元。對(duì)于模式物理量，其全球積分表示為：

(1)

因物理量在模式網(wǎng)格上離散化分布，設(shè)源網(wǎng)格中第個(gè)格點(diǎn)的量值為，面積為，上式可近似表示為：

(2)

設(shè)物理量的重映射為，則目標(biāo)網(wǎng)格與源網(wǎng)格滿足關(guān)系：

(3)

(4)

若該目標(biāo)網(wǎng)格共覆蓋源網(wǎng)格中個(gè)格點(diǎn)，則上式進(jìn)一步表示為：

(5)

其中：物理量通過下式進(jìn)行估計(jì)，

(6)

由于采用C網(wǎng)格離散方案,MPAS-A的水平動(dòng)量場(chǎng)定義在網(wǎng)格邊界并與其正交(Skamarock et al.,2012)。同時(shí),MPAS-A框架中采用徑向基函數(shù)(Bonaventura et al.,2011)將網(wǎng)格邊界的水平動(dòng)量轉(zhuǎn)換為網(wǎng)格中心的經(jīng)向風(fēng)和緯向風(fēng)分量,以滿足垂直動(dòng)量方程與參數(shù)化方案的需要。轉(zhuǎn)換接口在進(jìn)行風(fēng)場(chǎng)轉(zhuǎn)換時(shí),同樣使用的是網(wǎng)格中心的風(fēng)速分量。但值得注意的是,將GSI分析結(jié)果轉(zhuǎn)換至MPAS-A非結(jié)構(gòu)網(wǎng)格時(shí),必須將網(wǎng)格中心的風(fēng)速增量再投影至網(wǎng)格邊界,才能完成水平動(dòng)量場(chǎng)的更新。這一處理方案也與Ha et al.(2017)研究中所采用的方法一致。

1.2 垂直層坐標(biāo)轉(zhuǎn)換

在垂直方向,MPAS-A采用了跟隨地形的高度坐標(biāo)(Klemp,2011):

=+()(,,)。

(7)

其各個(gè)垂直層的氣壓通過狀態(tài)方程診斷獲得(Skamarock et al.,2012)。而GSI采用與業(yè)務(wù)模式GFS一致的地形追隨氣壓坐標(biāo)(Sela,2009),在其原始方案中,模式的三維氣壓場(chǎng)直接通過地表氣壓及兩組垂直廓線系數(shù)進(jìn)行確定:

(8)

因此,必須進(jìn)行必要的垂直坐標(biāo)轉(zhuǎn)換才能使得GSI系統(tǒng)內(nèi)部各個(gè)模式層上的氣象數(shù)據(jù)能夠正確對(duì)應(yīng)MPAS-A模式層上的氣壓和高度信息。

本研究采用反距離權(quán)重插值和線性外推方法,將MPAS-A垂直方向位于模式面半層的氣壓數(shù)據(jù)插值到整層,并輸入到GSI系統(tǒng)中進(jìn)行計(jì)算。這一方法避免了GSI系統(tǒng)采用其原始方案推導(dǎo)模式三維氣壓場(chǎng)。圖3為不同垂直層方案對(duì)于水汽混合比的處理結(jié)果與再分析資料的對(duì)比,其中圖3a為FNL分析場(chǎng),圖3b、圖3c分別表示基于GSI原始方案及本研究方案將圖3a所示變量處理至模式層后的垂直分布。從圖中可見,相比GSI原始方案(圖3b),本研究所用方案的處理結(jié)果在垂直方向與源數(shù)據(jù)保持一致(圖3c)。與氣壓類似,各模式層的高度信息也通過轉(zhuǎn)換接口向GSI系統(tǒng)傳遞。

圖3 不同垂直層方案下沿115°E的水汽混合比的垂直分布(單位:kg·kg-1):(a)FNL分析場(chǎng);(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 網(wǎng)格設(shè)置及轉(zhuǎn)換功能檢驗(yàn)

2.1 網(wǎng)格設(shè)置

本研究采用MPAS-A全球準(zhǔn)均勻網(wǎng)格,分辨率約30 km。根據(jù)譜截?cái)嗖〝?shù)的不同,GSI系統(tǒng)中預(yù)置了多套不同分辨率的高斯結(jié)構(gòu)網(wǎng)格。本文從中選取了T382方案,其緯向及經(jīng)向格點(diǎn)數(shù)分別為1 152和576,分辨率約35 km,與MPAS-A所用網(wǎng)格最為接近。在垂直方向,MPAS-A與GSI均設(shè)置為64層,其中MPAS-A的層頂高度為30 km。

2.2 網(wǎng)格轉(zhuǎn)換功能檢驗(yàn)

2.2.1 檢驗(yàn)方法

在GSI/MPAS同化及預(yù)報(bào)框架進(jìn)行應(yīng)用前,需首先對(duì)轉(zhuǎn)換接口進(jìn)行檢驗(yàn)。檢驗(yàn)時(shí)選取了2016年6月19日06時(shí)(世界時(shí),下同)的GFS 0.25°分析場(chǎng)產(chǎn)品,經(jīng)MPAS-A預(yù)處理模塊插值到非結(jié)構(gòu)網(wǎng)格作為基準(zhǔn)值,然后通過轉(zhuǎn)換接口將其轉(zhuǎn)換至GSI結(jié)構(gòu)網(wǎng)格后再轉(zhuǎn)換回原網(wǎng)格,通過以下誤差范數(shù)(Ullrich et al.,2009)衡量轉(zhuǎn)換結(jié)果同基準(zhǔn)值的差異:

(9)

(10)

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

(11)

和分別表示基準(zhǔn)值及經(jīng)過兩次轉(zhuǎn)換后的計(jì)算結(jié)果。表征全場(chǎng)大尺度的誤差特征，表征全場(chǎng)小尺度的誤差特征。

2.2.2 結(jié)果對(duì)比

圖4 模式物理量經(jīng)過兩次守恒重映射轉(zhuǎn)換后的誤差范數(shù)的垂直分布(a,b,c為L(zhǎng)1范數(shù);d,e,f為L(zhǎng)2范數(shù)):(a,d)緯向風(fēng);(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 模式物理量原始場(chǎng)(a,b;緯向風(fēng)單位:m/s;位溫單位:K)以及二階與一階精度守恒重映射轉(zhuǎn)換結(jié)果同原始場(chǎng)的絕對(duì)誤差百分率之間的差異(c,d;單位:%):(a,c)模式第25層緯向風(fēng);(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

在應(yīng)用一階或二階精度的守恒重映射方法進(jìn)行轉(zhuǎn)換時(shí),都需要利用預(yù)先生成的轉(zhuǎn)換權(quán)重系數(shù)進(jìn)行計(jì)算。在一次轉(zhuǎn)換過程中,采用一階精度方法由MPAS-A非結(jié)構(gòu)網(wǎng)格轉(zhuǎn)換至GSI結(jié)構(gòu)網(wǎng)格用時(shí)為93.996 18 s,從GSI結(jié)構(gòu)網(wǎng)格轉(zhuǎn)換回MPAS-A非結(jié)構(gòu)網(wǎng)格用時(shí)為47.220 25 s;二階精度方法在這兩個(gè)過程中的用時(shí)分別為119.549 10 s和65.299 86 s。從計(jì)算時(shí)效性方面來看,二階精度方法由于考慮了更多信息,在進(jìn)行轉(zhuǎn)換時(shí)相比一階方法消耗更多時(shí)間,但差異并不顯著。二階精度方法的耗時(shí)主要集中在轉(zhuǎn)換權(quán)重系數(shù)的生成環(huán)節(jié),在網(wǎng)格確定的情況下,只需預(yù)先生成轉(zhuǎn)換權(quán)重系數(shù)后便可重復(fù)使用。因此,采用一階或二階方法都能夠在實(shí)際運(yùn)行時(shí)實(shí)現(xiàn)快速轉(zhuǎn)換。

測(cè)試分析表明,采用二階精度的守恒重映射方法使得各變量的轉(zhuǎn)換誤差在大尺度和小尺度特征上相比一階精度都有減小,故二階精度守恒重映射方法對(duì)于變量的估計(jì)更加準(zhǔn)確,且耗時(shí)相比一階精度方法并未顯著增加。因此,為盡量減小轉(zhuǎn)換誤差對(duì)結(jié)果的影響,本研究對(duì)于各變量均采用二階精度守恒重映射方法進(jìn)行網(wǎng)格結(jié)構(gòu)轉(zhuǎn)換。

3 循環(huán)同化及預(yù)報(bào)試驗(yàn)

為進(jìn)一步分析GSI-MPAS框架的同化及預(yù)報(bào)結(jié)果,本節(jié)繼續(xù)基于30 km準(zhǔn)均勻網(wǎng)格,進(jìn)行了為期一周的連續(xù)循環(huán)同化及預(yù)報(bào)試驗(yàn),并將結(jié)果與ERA5再分析資料進(jìn)行對(duì)比。

3.1 試驗(yàn)設(shè)計(jì)

試驗(yàn)選取的時(shí)段為2016年6月19日00時(shí)—2016年6月26日00時(shí),初始場(chǎng)采用分辨率為0.25°的NCEP FNL數(shù)據(jù)?？刂圃囼?yàn)于起始時(shí)刻冷啟動(dòng)并進(jìn)行預(yù)報(bào),稱作“30 km_CTRL”。循環(huán)同化試驗(yàn)于2016年6月19日06時(shí)采用基于初始時(shí)刻的6 h預(yù)報(bào)場(chǎng)作為背景場(chǎng),開始進(jìn)行時(shí)間間隔為6 h的連續(xù)循環(huán)同化及預(yù)報(bào),稱作“30 km_DA”。對(duì)于00時(shí)或12時(shí)同化生成的分析場(chǎng),通過MPAS-A進(jìn)行24 h預(yù)報(bào),其他時(shí)刻同化完成后只進(jìn)行6 h預(yù)報(bào),以獲得下一分析時(shí)刻的背景場(chǎng)。

MPAS-A模式為v7.0版本,采用Grell-Freitas積云對(duì)流方案(Grell and Freitas,2014),WSM6微物理過程方案(Hong and Lim,2006)、YSU邊界層方案(Hong et al.,2006)、RRTMG長(zhǎng)波及短波方案(Iacono et al.,2008)以及Noah陸面過程方案(Chen and Dudhia,2001)等。另外,在試驗(yàn)期間啟用了海溫更新,模式每逢00時(shí)從與初始場(chǎng)規(guī)格一致的FNL資料中獲取相應(yīng)時(shí)刻的海溫及海冰數(shù)據(jù)。

GSI系統(tǒng)采用全球三維變分分析方案,其控制變量為流函數(shù)、非平衡速度勢(shì)、非平衡虛溫、非平衡地表氣壓、標(biāo)準(zhǔn)化相對(duì)濕度以及臭氧和云凝結(jié)混合比(Kleist et al.,2009),采用官方提供的全球背景誤差協(xié)方差。觀測(cè)數(shù)據(jù)為GDAS全球常規(guī)觀測(cè)資料,其2016年6月19日12時(shí)的觀測(cè)種類及位置分布如圖6所示。

圖6 循環(huán)同化試驗(yàn)于2016年6月19日12時(shí)采用的觀測(cè)資料種類(括號(hào)內(nèi)數(shù)字為觀測(cè)數(shù)量)及其分布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 同化及預(yù)報(bào)結(jié)果評(píng)估

3.2.1 高空及地面分析場(chǎng)

各組試驗(yàn)在最后一個(gè)同化時(shí)刻(2016年6月26日00時(shí))的500 hPa位勢(shì)高度與溫度分析場(chǎng)與ERA5資料的差異如圖7所示。從圖中可見,控制試驗(yàn)?zāi)M的500 hPa溫度場(chǎng)相比ERA5資料在中緯度地區(qū)有明顯偏高,而同化試驗(yàn)在一定程度上減小了控制試驗(yàn)的溫度偏差,以北半球亞歐大陸南部最為顯著。圖7d為該同化時(shí)刻的500 hPa溫度分析增量場(chǎng),可見由于北半球觀測(cè)資料較多(圖6),分析增量主要集中在北半球地區(qū),使得該地區(qū)的改進(jìn)相比南半球更大。從高度場(chǎng)上看,循環(huán)同化試驗(yàn)的模擬結(jié)果相比控制試驗(yàn)更接近ERA5資料,特別是對(duì)北半球的槽脊系統(tǒng)定位更加準(zhǔn)確。圖7e進(jìn)一步描繪了循環(huán)同化試驗(yàn)及控制試驗(yàn)的北半球500 hPa高度場(chǎng)與ERA5資料的距平相關(guān)系數(shù)隨模擬時(shí)間的變化。從此圖可見,控制試驗(yàn)的高度場(chǎng)距平相關(guān)系數(shù)隨模擬時(shí)間的延伸而不斷降低,而循環(huán)同化試驗(yàn)由于不斷吸收觀測(cè)數(shù)據(jù),其高度場(chǎng)距平相關(guān)系數(shù)減小幅度明顯低于控制試驗(yàn),并在相當(dāng)長(zhǎng)的時(shí)間內(nèi)保持在較高水平。

圖7 2016年6月26日00時(shí)500 hPa溫度場(chǎng)(陰影;單位:K)和高度場(chǎng)(等值線;單位:gpm):(a)控制試驗(yàn);(b)循環(huán)同化試驗(yàn);(c)ERA5資料;(d)該時(shí)刻500 hPa溫度場(chǎng)的分析增量;(e)7 d試驗(yàn)期間500 hPa高度場(chǎng)與ERA5資料的距平相關(guān)系數(shù)隨時(shí)間的變化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為各組試驗(yàn)的平均海平面氣壓模擬結(jié)果與ERA5資料的差異。可以看出,控制試驗(yàn)結(jié)果中存在的偏差大值區(qū)較多(圖8a),反映出其對(duì)于主要地面天氣系統(tǒng)的模擬不夠準(zhǔn)確。而在同化試驗(yàn)中,偏差水平在北半球有明顯降低,整體體現(xiàn)為微弱的正偏差,在南半球依然存在少量偏差大值區(qū),可能仍與南半球同化的觀測(cè)資料數(shù)量偏少有關(guān)。

圖8 2016年6月26日00時(shí)平均海平面氣壓場(chǎng)模擬結(jié)果與ERA5資料的差異(單位:hPa):(a)控制試驗(yàn);(b)循環(huán)同化試驗(yàn)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

從整體形勢(shì)場(chǎng)上看,GSI-MPAS同化及預(yù)報(bào)框架能夠有效吸收多源觀測(cè)信息,并對(duì)溫壓場(chǎng)起到合理的調(diào)整作用。

3.2.2 同化及預(yù)報(bào)結(jié)果的均方根誤差

為定量分析GSI-MPAS同化及預(yù)報(bào)框架的分析場(chǎng)及預(yù)報(bào)場(chǎng)特征,將循環(huán)同化及預(yù)報(bào)試驗(yàn)在各個(gè)00時(shí)與12時(shí)的溫度和緯向風(fēng)分析場(chǎng)及24 h預(yù)報(bào)場(chǎng)分別與ERA5資料在不同區(qū)域計(jì)算均方根誤差并進(jìn)行時(shí)間平均(圖9)。從全球整體結(jié)果上看,同化試驗(yàn)的均方根誤差在各個(gè)氣壓層相比控制試驗(yàn)都有減小,對(duì)于控制試驗(yàn)誤差較大的氣壓層,改進(jìn)作用也更加明顯。從不同地區(qū)結(jié)果上看,兩組試驗(yàn)在北半球的模擬誤差整體低于南半球,特別是大氣中低層;而且同化試驗(yàn)相比控制試驗(yàn)在北半球帶來的改進(jìn)幅度也大于南半球。在熱帶地區(qū),控制試驗(yàn)在對(duì)流層中高層預(yù)報(bào)誤差較大,但由于觀測(cè)資料稀少,同化試驗(yàn)在該地區(qū)的改進(jìn)作用也最小。24 h預(yù)報(bào)場(chǎng)中,各變量模擬誤差的整體分布特征與分析場(chǎng)比較一致,但預(yù)報(bào)場(chǎng)的誤差量值相比分析場(chǎng)有所增加。同時(shí),同化試驗(yàn)在南北半球的模擬誤差差異也有進(jìn)一步擴(kuò)大。兩組試驗(yàn)的均方根誤差評(píng)分結(jié)果進(jìn)一步表明,高空觀測(cè)資料的引入對(duì)于模式三維變量場(chǎng)具有整體改進(jìn)作用,在觀測(cè)資料分布密集的北半球地區(qū)帶來的改進(jìn)更為顯著。

圖9 溫度(a,d;單位:K)、緯向風(fēng)(b,e;單位:m/s)和相對(duì)濕度(c,f;單位:%)的模擬結(jié)果在不同地區(qū)相比于ERA5資料的均方根誤差廓線:(a,b,c)分析場(chǎng);(d,e,f)24 h預(yù)報(bào)場(chǎng)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 我國(guó)東部梅雨期降水模擬

因循環(huán)同化試驗(yàn)期間正值2016年梅雨期(趙俊虎等,2018),圖10比較了同化預(yù)報(bào)時(shí)段內(nèi)華東地區(qū)(114°～123°E,28°～36°N)梅雨降水的緯向平均值隨時(shí)間的變化,這能夠反映出主要雨帶的南北移動(dòng)情況。其中,降水實(shí)況為自動(dòng)站與CMORPH降水產(chǎn)品融合的逐時(shí)降水?dāng)?shù)據(jù)集。從圖中可見,控制試驗(yàn)在初期模擬的雨帶相比觀測(cè)有明顯北移,而同化試驗(yàn)在一定程度上消除了控制試驗(yàn)的過度空?qǐng)?bào),盡管其降水強(qiáng)度相比觀測(cè)有所減小。隨著模擬時(shí)間的推移,控制試驗(yàn)對(duì)于降水的預(yù)測(cè)能力出現(xiàn)明顯降低,其模擬的降水中心與實(shí)況偏離甚遠(yuǎn)。同化試驗(yàn)仍然能夠捕捉到6月21日及23日出現(xiàn)的強(qiáng)降水事件,雖然在降水位置及強(qiáng)度上與實(shí)況仍有差距,但相比控制試驗(yàn)改進(jìn)明顯,顯示出同化觀測(cè)資料對(duì)于改善降水預(yù)報(bào)的積極影響。

圖10 中國(guó)東部梅雨帶(114°～123°E,28°～36°N)的緯向平均降水隨時(shí)間的變化(單位:mm/h):(a)實(shí)況;(b)控制試驗(yàn);(c)循環(huán)同化試驗(yàn)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 結(jié)論與討論

為改進(jìn)基于非結(jié)構(gòu)計(jì)算網(wǎng)格的全球非靜力大氣模式MPAS-A的模擬性能,本研究通過基于守恒重映射方法的球面網(wǎng)格轉(zhuǎn)換和垂直坐標(biāo)轉(zhuǎn)換,引入GSI全球三維變分同化系統(tǒng),形成了GSI-MPAS同化及預(yù)報(bào)框架。隨后,開展了一系列數(shù)值試驗(yàn)及循環(huán)同化預(yù)報(bào)試驗(yàn),得到如下結(jié)論:

1)守恒重映射方法能夠?qū)崿F(xiàn)氣象要素在MPAS-A的非結(jié)構(gòu)網(wǎng)格與GSI的高斯結(jié)構(gòu)網(wǎng)格之間的轉(zhuǎn)換。二階精度的守恒重映射方法相比一階精度方法的轉(zhuǎn)換誤差更小,能對(duì)物理量場(chǎng)進(jìn)行更加準(zhǔn)確的轉(zhuǎn)換。

2)經(jīng)過GSI同化后,MPAS-A預(yù)報(bào)得到的動(dòng)力場(chǎng)和熱力場(chǎng)結(jié)果在整體上獲得改進(jìn),表明多源觀測(cè)資料在GSI-MPAS同化及預(yù)報(bào)框架中得到了有效應(yīng)用。

3)對(duì)試驗(yàn)結(jié)果的初步分析表明,在觀測(cè)資料分布較為密集的北半球地區(qū),同化所得分析場(chǎng)帶來的改進(jìn)作用更加明顯,其各個(gè)變量的預(yù)報(bào)結(jié)果也更為合理。同時(shí),華東地區(qū)梅雨期的降雨預(yù)報(bào)對(duì)比也表明同化后的降水預(yù)報(bào)結(jié)果在時(shí)間和空間上均與觀測(cè)更加接近。

從目前初步應(yīng)用結(jié)果可見,本研究構(gòu)建的GSI-MPAS同化及預(yù)報(bào)框架能夠利用多源觀測(cè)資料得到更加合理的初始場(chǎng)。但通過比較不同地區(qū)的同化及預(yù)報(bào)結(jié)果也表明,因本文同化的常規(guī)觀測(cè)資料大多

數(shù)位于北半球,導(dǎo)致南半球及熱帶地區(qū)的同化改進(jìn)相對(duì)較小。因此,下一步將考慮同化更多觀測(cè),如衛(wèi)星資料(陳耀登等,2021)等,以減小南北半球分析結(jié)果的差異,使全球整體分析結(jié)果更加合理。另外,GSI系統(tǒng)的背景誤差及相關(guān)資料的觀測(cè)誤差參數(shù)設(shè)置也需進(jìn)行進(jìn)一步優(yōu)化。

本論文的數(shù)值計(jì)算得到了南京信息工程大學(xué)高性能計(jì)算中心的計(jì)算支持和幫助。