程銳 宇如聰 徐幼平 劉娟 黃靜

1 地理信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安710054

2 中國(guó)氣象局，北京 100081

3 中國(guó)科學(xué)院大氣物理研究所大氣科學(xué)和地球流體力學(xué)數(shù)值模擬國(guó)家重點(diǎn)實(shí)驗(yàn)室(LASG)，北京 100029

1 引言

開(kāi)發(fā)動(dòng)力框架和設(shè)計(jì)相應(yīng)的求解算法是模式發(fā)展的焦點(diǎn)問(wèn)題之一。WRF（Weather Research and Forecasting model）動(dòng)力框架考慮了更高階數(shù)值精度和標(biāo)量守恒為特性的數(shù)值算法，物理過(guò)程選項(xiàng)豐富、計(jì)算高效，同時(shí)具備可移植性、并行化、可擴(kuò)展性；用于封裝模式的軟件框架考慮了多層級(jí)軟件架構(gòu)以及應(yīng)用程序編程接口（API）設(shè)計(jì)等功能（Michalakes et al., 2001），可組裝動(dòng)力框架、物理過(guò)程、I/O 等模式元素?？梢钥闯?，該動(dòng)力框架具備設(shè)計(jì)先進(jìn)性、算法精確性和求解高效性，既能把模式元素富有邏輯、高效科學(xué)地組織起來(lái)，同時(shí)也確保了系統(tǒng)的可擴(kuò)展、易于開(kāi)發(fā)及大規(guī)模并行計(jì)算等功能，是現(xiàn)今利用率較高的先進(jìn)模式框架。

WRF 模式在東亞地區(qū)獲得了廣泛應(yīng)用，為中尺度天氣預(yù)報(bào)機(jī)理研究和暴雨、臺(tái)風(fēng)等模擬預(yù)報(bào)提供了有力工具；但和其他數(shù)值模式一樣，對(duì)于復(fù)雜地形下暴雨預(yù)報(bào)性能仍然不高。借鑒Eta 坐標(biāo)模式對(duì)陡峭地形有效處理（Mesinger, 1984, 2004; 宇如聰 和 徐 幼 平, 2004; 普 業(yè) 等, 2008; 程 銳 等, 2018,2019），若能在WRF 模式中采用階梯地形垂直坐標(biāo)進(jìn)行復(fù)雜地形表達(dá)，改善陡峭地形下氣壓梯度力、水平平流和擴(kuò)散等過(guò)程的刻畫(huà)，則既是對(duì)該模式功能的有益提升，也能夠?qū)鉀Q復(fù)雜地形下暴雨預(yù)報(bào)問(wèn)題提供一定參考。

1.1 垂直坐標(biāo)設(shè)計(jì)

垂直坐標(biāo)構(gòu)造是數(shù)值模式的重要內(nèi)容；理論上，在數(shù)值天氣預(yù)報(bào)中，凡與幾何高度有單調(diào)關(guān)系的任何變量都可作為垂直坐標(biāo)，不同坐標(biāo)變換僅僅帶來(lái)方程組形式變化以及邊界條件的改變。Richardson最初進(jìn)行數(shù)值天氣預(yù)報(bào)設(shè)計(jì)時(shí)采用了幾何高度坐標(biāo)，并由Kasahara and Washington（1967）實(shí)現(xiàn)利用高度坐標(biāo)系進(jìn)行數(shù)值天氣預(yù)報(bào)。Eliassen（1949）提出氣壓作為垂直坐標(biāo)的概念，并得到廣泛應(yīng)用；后來(lái)，為了更好處理下邊界，Phillips（1957）使用了歸一化坐標(biāo)（σ=p/ps）；再后來(lái)，由歸一化坐標(biāo)逐漸發(fā)展為基于高度（Gal-Chen and Somerville, 1975）和質(zhì)量（Laprise, 1992）的地形追隨坐標(biāo)?；贓gger（1972）的地形障礙方法，曾慶存1982 年提出一種改進(jìn)的Sigma 坐標(biāo)（手稿；稿件題目：修改的σ 坐標(biāo)），與Mesinger（1984）提出的Eta 坐標(biāo)一致，都為階梯地形坐標(biāo)?，F(xiàn)在比較流行的垂直混合坐標(biāo)的概念由Sangster（1960）提出，他建議模式低層采用Sigma 坐標(biāo)，高層逐漸采用氣壓坐標(biāo)，從而在高層避免氣壓梯度力大量小差問(wèn)題。Steppeler et al.（2002）將模式地形表示為線性連續(xù)樣條函數(shù)形式，相當(dāng)于傾斜階梯地形坐標(biāo)。

在選擇垂直坐標(biāo)時(shí)，首先要考慮是否適合處理地形。地形是大氣運(yùn)動(dòng)重要的強(qiáng)迫源，既可將其看作阻滯氣流的巨大障礙物，也可看作是矗立在對(duì)流層的熱島和水汽供應(yīng)源（紀(jì)立人等, 2005）；地形還是對(duì)流活動(dòng)的觸發(fā)源。中尺度數(shù)值模式的主要目的是描寫(xiě)大氣中與災(zāi)害天氣有關(guān)的強(qiáng)對(duì)流活動(dòng)，而大陸上的對(duì)流災(zāi)害天氣又往往和地形強(qiáng)迫作用密切相關(guān)。另外，隨著模式分辨率的提高，陡峭地形坡度越來(lái)越大；比如2 km 分辨率的中尺度模式，地形坡度能達(dá)到0.3 左右，而一般的天氣模式地形坡度小于0.01（胡江林和王盤(pán)興, 2007）。因此，精確、合理的構(gòu)造模式地形是數(shù)值天氣預(yù)報(bào)的重要問(wèn)題，尤其我國(guó)境內(nèi)存在青藏高原等復(fù)雜地形，與地形相關(guān)的數(shù)值計(jì)算與處理更值得深入研究。

1.2 地形追隨坐標(biāo)的問(wèn)題及其改進(jìn)

地形追隨坐標(biāo)下，因坐標(biāo)系非正交，所以會(huì)產(chǎn)生虛假曲率影響；由于坐標(biāo)面不水平，水平平流和擴(kuò)散處理需要額外人為控制；尤其是陡峭地形下氣壓梯度力的計(jì)算精度不高（Mesinger, 2004）。當(dāng)然，地形追隨坐標(biāo)對(duì)陡峭地形處氣壓梯度力的求解誤差可帶來(lái)氣象場(chǎng)模擬預(yù)報(bào)偏差（李興良等, 2005;李興良和陳德輝, 2015）。選用該垂直坐標(biāo)，對(duì)背風(fēng)坡氣旋、冷空氣阻塞等山脈下風(fēng)方天氣事件的模擬和預(yù)報(bào)也存在困難。

對(duì)于地形追隨坐標(biāo)模式，提高水平氣壓梯度力計(jì)算精度有多種應(yīng)對(duì)方案，包括設(shè)計(jì)協(xié)調(diào)的差分格式（Corby et al., 1972；該格式在MM5 中得到應(yīng)用）、反插回等壓面計(jì)算（Smagorinsky et al., 1967；胡江林和王盤(pán)興, 2007）、扣除法（曾慶存，1963，1979；錢(qián)永甫和周天軍, 1995；WRF 模式中也采用了這種方法）以及上下層平均法（錢(qián)永甫等, 1978）等。另外，可將σ坐標(biāo)變換為p–σ混合坐標(biāo)，從而提高氣象要素模擬精度。事實(shí)上，錢(qián)永甫（1985）就提出一個(gè)5 層p–σ混合坐標(biāo)原始方程模式；趙鳴和方娟（1996）將其發(fā)展為9 層混合坐標(biāo)模式。Steppeler et al.（2002）將模式地形表示為線性連續(xù)樣條函數(shù)形式，并使用z和σ混合坐標(biāo)構(gòu)造LMz（z-coordinate nonhydrostatic version of the Lokal model）模式，其中動(dòng)力框架在z坐標(biāo)計(jì)算而物理過(guò)程在σ坐標(biāo)計(jì)算。該方法是地形追隨坐標(biāo)和階梯地形坐標(biāo)的折衷方案，可以同時(shí)克服上述階梯地形坐標(biāo)及追隨地形坐標(biāo)的一些問(wèn)題。WRF 模式也將混合坐標(biāo)開(kāi)發(fā)作為模式更新的重要內(nèi)容，在其混合垂直坐標(biāo)框架中，靠近地形的模式低層采用地形追隨坐標(biāo)，而地形之上逐漸向等壓面坐標(biāo)過(guò)渡松弛；這樣既滿足了易于地形處理的需要，同時(shí)也減小了地形對(duì)于模式高層的人為影響，從而使對(duì)流層頂?shù)哪Ｊ酱硇缘玫皆鰪?qiáng)。

上述修正主要針對(duì)動(dòng)力計(jì)算層面，還可以通過(guò)減緩模式坐標(biāo)面坡度，從而減小地形追隨坐標(biāo)計(jì)算誤差的方式提高計(jì)算精度。Sch?r et al.（2002）提出一種新的平滑層垂直坐標(biāo)（SLEVE），他們通過(guò)平滑復(fù)雜地形以上Sigma 坐標(biāo)面來(lái)提高計(jì)算精度，并成功將此方法運(yùn)用于COSMO（Consortium for Small-scale Modeling）模式中；后來(lái)，Z?ngl（2003）將該方法推廣到質(zhì)量追隨坐標(biāo)模式中。美國(guó)跨尺度預(yù)報(bào)模式（Model for Prediction Across Scales，簡(jiǎn)稱(chēng)MPAS）采用平滑地形追隨（STF）混合坐標(biāo)（Klemp, 2011），通過(guò)控制地形對(duì)坐標(biāo)面的影響使得高層坐標(biāo)面逐漸水平，從而使氣流更趨平滑，這種調(diào)制作用在模式分辨率提高時(shí)尤顯重要。相比SLEVE 方法，STF 方法對(duì)地形追隨坐標(biāo)面平滑更大。李超等（2019）基于GRAPES 模式對(duì)國(guó)際先進(jìn)的平緩-混合坐標(biāo)（Smoothed and Hybrid Terrain Following coordinate）進(jìn)行了系統(tǒng)分析、試驗(yàn)，設(shè)計(jì)一種改進(jìn)的余弦三角函數(shù)為基函數(shù)的平緩-混合坐標(biāo)（COS 坐標(biāo)），使低層坐標(biāo)面垂直分布更加均勻，低層地形作用衰減的垂直變化更加均勻，減小了計(jì)算誤差，提高了計(jì)算穩(wěn)定性。張旭等人（2015）也利用GRAPES 區(qū)域模式，開(kāi)展了Klemp（2011）和Sch?r et al.（2002）高度地形追隨坐標(biāo)理想試驗(yàn)研究，他們主要通過(guò)設(shè)定不同的地形衰減廓線或?qū)ψ鴺?biāo)面地形高度進(jìn)行平滑來(lái)減小小尺度地形影響。屠妮妮等（2012）、何光碧等（2015）通過(guò)濾波方案直接對(duì)模式地形進(jìn)行處理，減緩地形坡度，提高模式預(yù)報(bào)性能。

為了解決追隨地形坐標(biāo)的非正交問(wèn)題，很多人都嘗試在陡峭地形之上使用數(shù)值網(wǎng)格生成方法設(shè)計(jì)正交的地形追隨網(wǎng)格，從而減小模式截?cái)嗾`差，如Erdun（1997）提出Schwarz-Christoffel 正形變換方法并在復(fù)雜地形上產(chǎn)生2D 正交網(wǎng)格，Li et al.（2014）則設(shè)計(jì)提出3D 正交曲線地形追隨坐標(biāo)。

隨著計(jì)算能力大幅提升，Eta 坐標(biāo)模式在邊界層內(nèi)分層也能得到精細(xì)表述，對(duì)重力波和邊界層特征的刻畫(huà)能力得到提高，加之其對(duì)水平氣壓梯度力和平流的刻畫(huà)優(yōu)勢(shì)，從而將Eta 坐標(biāo)的應(yīng)用重新提上模式發(fā)展日程。Mesinger et al.（2012）設(shè)計(jì)了一種傾斜Eta 坐標(biāo)，較好清除了Gallus and Klemp（2000）發(fā)現(xiàn)的在原Eta 階梯地形模式中存在的階梯轉(zhuǎn)角的虛假渦度。隨后，Mesinger and Veljovic（2017）進(jìn)行了該版Eta 坐標(biāo)模式和Sigma 坐標(biāo)模式比較試驗(yàn)，發(fā)現(xiàn)當(dāng)Eta 垂直分層加密后氣流分離現(xiàn)象得到消除。另外，Eta 模式在云物理和輻射過(guò)程以及資料同化改進(jìn)后，其在降水評(píng)分上相對(duì)于Sigma 模式保持較高水平；Eta 模式優(yōu)勢(shì)在美國(guó)西部特別明顯，對(duì)于所有降水類(lèi)型Eta 模式性能更加優(yōu)異，特別是較強(qiáng)降水。

1.3 山脈波的求解

山脈波的經(jīng)典求解是給定穩(wěn)定層結(jié)和平直西風(fēng)情形下，分析鐘形山脈激發(fā)的重力波解。實(shí)際鐘形山脈可以進(jìn)行階梯地形近似，這種近似的準(zhǔn)確性主要與地形的水平尺度以及垂直分辨率有關(guān)。1943年Lyra 推導(dǎo)了帽型山脈氣流過(guò)山的分析解（Alaka,1960）。Gallus et al.（2000）在此基礎(chǔ)上重新推導(dǎo)了旋轉(zhuǎn)大氣的山脈波解。

對(duì)于靜力效應(yīng)顯著的波動(dòng)，波動(dòng)能量和動(dòng)量垂直傳播，并且波動(dòng)局限于山脈上方區(qū)域，定常位相隨高度向上游傾斜。此時(shí)，在轉(zhuǎn)角之上產(chǎn)生的擾動(dòng)的尺度比山脈強(qiáng)迫的尺度小好多，因此會(huì)產(chǎn)生明顯的歪曲現(xiàn)象，也就使得每個(gè)階梯轉(zhuǎn)角之上的環(huán)流表現(xiàn)出強(qiáng)局地特征。雖然轉(zhuǎn)角之上的擾動(dòng)振幅與山脈寬度無(wú)關(guān)，但垂直速度隨著山脈寬度增加而減小。因此，當(dāng)山脈寬度增加時(shí)，這些階梯轉(zhuǎn)角擾動(dòng)的相對(duì)振幅變得愈加明顯。

2 Eta 坐標(biāo)下WRF 動(dòng)力框架

從原始方程出發(fā)，推導(dǎo)考慮地圖投影變換和擾動(dòng)變換的Eta 坐標(biāo)WRF 濕動(dòng)力框架如下：

（a）第一運(yùn)動(dòng)方程：

（b）第二運(yùn)動(dòng)方程：

（c）第三運(yùn)動(dòng)方程：

（d）連續(xù)方程：

（e）熱力學(xué)方程：

（f）擾動(dòng)位勢(shì)方程：

（g）水物質(zhì)方程：

（i）靜力近似方程：

上述9 個(gè)方程共同組成了Eta 坐標(biāo)下WRF 動(dòng)力框架。其中，上標(biāo)“ ′”表示擾動(dòng)量，下標(biāo)“d”、“m”分別表示干、濕過(guò)程； γ=cp/cv，cp、cv和Rd分別表示定壓比熱、定容比熱和干空氣氣體常數(shù)；qv+qc+qr+qi+···表示水汽和水凝物混合比之和；(U,V,W)=(Pu,Pv,Pw)，?=P，Θ=Pθ，Qm=Pqm，qm表 示水凝物混合比；FU、FV、FW主要為科氏力、曲率項(xiàng)力、混合項(xiàng)和物理過(guò)程轉(zhuǎn)換項(xiàng)的貢獻(xiàn)；FΘ、FQm主要為混合項(xiàng)和物理過(guò)程轉(zhuǎn)換項(xiàng)的貢獻(xiàn)；η表示Eta 坐標(biāo)：

其中，

式中， π表示靜力氣壓， πt為模式頂靜力氣壓， πs為地面靜力氣壓， ηs表 示地形處Eta 值； πrf為參考靜力氣壓；zs為 地形高度，zb為 平緩地形高度， α表示大氣比容，p0表示參考?xì)鈮海话闳?000 hPa。

其他變量使用氣象常用符號(hào)表示，不再贅述。

3 Eta 垂直坐標(biāo)設(shè)計(jì)

3.1 參考狀態(tài)定義

標(biāo)準(zhǔn)層結(jié)近似由曾慶存先生提出，宇如聰（1989）在中尺度模式中應(yīng)用該近似方法實(shí)現(xiàn)靜力扣除計(jì)算。本文利用標(biāo)準(zhǔn)層結(jié)大氣構(gòu)造階梯地形，具體形式如下：

3.2 階梯地形表征

3.2.1 階梯地形2D 表征量

由于階梯地形坐標(biāo)具有準(zhǔn)水平的性質(zhì)，因此坐標(biāo)面和地形存在交割關(guān)系，在不同的水平格點(diǎn)，大氣變量在垂直方向的分布將會(huì)有三種情況：地形之下、地形表面或地形之上。在設(shè)計(jì)模式積分計(jì)算時(shí)，地形以下層次將不能參與。考慮WRF 水平網(wǎng)格分布為Arakawa-C 型跳點(diǎn)的特點(diǎn)，設(shè)計(jì)變量ksh（非界面層質(zhì)量點(diǎn)緊鄰地形的垂直層次編號(hào)）、ksu（緯向風(fēng)點(diǎn)緊鄰地形的垂直層次編號(hào)）和ksv（經(jīng)向風(fēng)點(diǎn)緊鄰地形的垂直層次編號(hào)），對(duì)緊貼地形的Eta 坐標(biāo)層次進(jìn)行判斷。模式垂直方向采用Lorenz 跳點(diǎn)形式，界面層和中間層交替分布，垂直速度和位勢(shì)高度分布在界面層，其它變量分布在中間層，且從地表至模式層頂，垂直層次編號(hào)逐漸遞增。

圖1（a）Eta-WRF 模式的階梯地形高度（單位：m）及階梯坐標(biāo)模式2D 表征量（b）ksh、（c）ksu、（d）ksv 隨水平網(wǎng)格的變化Fig.1 (a) Eta-WRF Stepped mountain height (units: m) and variation of Eta model 2D-mesh indicators with horizontal grid: (b) ksh, (c) ksu, and (d)ksv

下面，我們將通過(guò)簡(jiǎn)單階梯地形來(lái)分析ksh、ksu 和ksv 等2D 表征量與階梯地形之間的關(guān)系（圖1）。在試驗(yàn)中，鐘形山脈實(shí)際最高高度為100 m，山脈半寬10 km，尺度系數(shù)選為zscale=0.5。從圖中可以看到，在區(qū)域中心存在一個(gè)孤立地形，高約100 m；與階梯地形對(duì)應(yīng)的區(qū)域中心附近ksh、ksu 和ksv 都為2，而在其它區(qū)域?yàn)?。我們以質(zhì)量點(diǎn)（不包括界面層變量垂直速度和位勢(shì)高度）表征量ksh 為例來(lái)說(shuō)明：通過(guò)階梯地形設(shè)計(jì)，可使階梯地形表面正好處于界面層，且緊鄰階梯地形之上的中間層和該界面層具有相同垂直層次編號(hào)；因此，“ksh=2”表示區(qū)域中心附近質(zhì)量點(diǎn)階梯地形所在層次為模式第2 個(gè)界面層。

3.2.2 3D 表征量

有了上述階梯地形2D 表征量，就能夠開(kāi)展階梯地形坐標(biāo)下動(dòng)力計(jì)算。但對(duì)于3D 格點(diǎn)空間的計(jì)算問(wèn)題，首先進(jìn)行水平網(wǎng)格判斷，再進(jìn)行垂直方向計(jì)算，會(huì)破壞數(shù)組的存取規(guī)則，影響計(jì)算效率。因此，為了方便高效開(kāi)展Eta 垂直坐標(biāo)模式3D 網(wǎng)格計(jì)算，需引入階梯地形3D 表征量mph、mphi、mpu、mpv，分別表示模式中間層質(zhì)量點(diǎn)所在3D 網(wǎng)格與階梯地形關(guān)系、界面層質(zhì)量點(diǎn)所在3D 網(wǎng)格與階梯地形關(guān)系、緯向速度點(diǎn)所在3D 網(wǎng)格與階梯地形關(guān)系以及經(jīng)向速度點(diǎn)所在3D 網(wǎng)格與階梯地形關(guān)系。

圖2 的計(jì)算方案與圖1 相同，為了更方便刻畫(huà)階梯地形3D 表征變量與Eta 垂直層次的關(guān)系，以及3D 表征量和2D 表征量表達(dá)的一致性，我們利用圖中填充的陰影部分表示階梯地形以下的網(wǎng)格點(diǎn)（其值設(shè)置為0）；此處Y軸表示Eta 層次，而不是3D 表征量的具體數(shù)值。從圖中可以發(fā)現(xiàn)，在區(qū)域中心階梯地形處，半層質(zhì)量點(diǎn)和速度點(diǎn)第2 層以下mph、mpu、mpv 為0，整層質(zhì)量點(diǎn)第3 層以下mphi 為0。再以mph 為例進(jìn)行說(shuō)明：mph 在中心區(qū)域以外沒(méi)有填充，表示中心區(qū)域以外沒(méi)有階梯地形存在，這和圖1a 是相互對(duì)應(yīng)的；在中心附近區(qū)域，因?yàn)榇嬖诠铝㈦A梯地形，且其位于第2 個(gè)界面層（ksh=2），因此中心區(qū)域附近第2 個(gè)界面層以下mph=0。對(duì)于界面層表征量mphi（圖2b），因?yàn)橹行母浇鼌^(qū)域階梯地形位于第2 個(gè)界面層，因此第3 個(gè)界面層以上才能保證存在階梯地形以上網(wǎng)格點(diǎn)，也就是說(shuō)第3 個(gè)界面層以下mphi=0；而在中心以外區(qū)域無(wú)階梯地形，故第2 個(gè)界面層以上就能保證存在階梯地形以上網(wǎng)格點(diǎn)，也就是說(shuō)第2 個(gè)界面層以下mphi=0。

4 山脈波試驗(yàn)

4.1 試驗(yàn)方案

原WRF（即Sigma 坐標(biāo)下動(dòng)力框架）山脈波試驗(yàn)方案作為對(duì)照試驗(yàn)（見(jiàn)表1），而將垂直坐標(biāo)變換和積分時(shí)間延長(zhǎng)作為比對(duì)試驗(yàn)（具體方案見(jiàn)表2）。通過(guò)數(shù)值試驗(yàn)，一方面檢驗(yàn)更新的模式動(dòng)力框架的正確性和有效性；另一方面通過(guò)對(duì)比試驗(yàn)考查垂直坐標(biāo)選擇對(duì)地形波模擬的影響。我們假定山脈地形最高為100 m，鐘形山脈最高處位于區(qū)域中心，山脈半寬為10 個(gè)格距（20 km，圖略），Brunt-Vaisala 頻率設(shè)為10?2s?1。當(dāng)前，大氣背景場(chǎng)設(shè)置為穩(wěn)定層結(jié)干大氣，水平平直西風(fēng)10 m s?1。數(shù)值模式選用2D 模型，模式區(qū)域水平延展202 個(gè)格點(diǎn)，水平分辨率為2 km，垂直方向指數(shù)拉伸分80 層，時(shí)間步長(zhǎng)20 s，模擬時(shí)長(zhǎng)36000 s。X方向?yàn)檩椛溟_(kāi)邊界，Y方向?yàn)橹芷谶吔?，Z方向剛體邊界條件。水平平流采用5 階方案計(jì)算，垂直平流采用3 階方案計(jì)算。當(dāng)前試驗(yàn)不考慮任何物理過(guò)程。不難看出，本方案模擬對(duì)象為典型靜力特性的山脈波動(dòng)，對(duì)這類(lèi)波動(dòng)的數(shù)值模擬也是Eta 坐標(biāo)模式的難點(diǎn)問(wèn)題之一。

表1 原WRF 模式（Sigma 坐標(biāo)）山脈波試驗(yàn)方案（對(duì)照試驗(yàn)）Table1 Original WRF (Sigma coordinates) configuration for mountain wave simulation (control experiment)

表2 Eta 坐標(biāo)對(duì)比試驗(yàn)方案Table2 Configuration for comparative experiments with Eta coordinates

4.2 對(duì)照試驗(yàn)結(jié)果分析

圖3 給出Sigma 坐標(biāo)下WRF 動(dòng)力框架對(duì)山脈波的模擬（對(duì)照試驗(yàn)EXP01）。從圖1 可以發(fā)現(xiàn)，積分36000 s 時(shí)，由于鐘形山脈的動(dòng)力強(qiáng)迫，出現(xiàn)了明顯的地形駐波特征。從溫度擾動(dòng)可以發(fā)現(xiàn)，在鐘形山脈山脊靠近背風(fēng)坡一側(cè)從低層到高層存在3個(gè)大值中心（最大約0.45 K，位于模式第10 層左右），強(qiáng)度隨高度減弱。氣壓擾動(dòng)呈現(xiàn)高層弱而低層強(qiáng)、迎風(fēng)坡弱而背風(fēng)坡強(qiáng)的形勢(shì)，最強(qiáng)擾動(dòng)中心在背風(fēng)坡靠近山脊附近，并位于模式底層，強(qiáng)度約?12 Pa。垂直速度分布和擾動(dòng)位溫相似，在山脊靠近背風(fēng)坡一層從低層至高層存在交錯(cuò)分布的上升、下沉運(yùn)動(dòng)帶，最大強(qiáng)度可達(dá)到0.04 m s?1，中心強(qiáng)度隨高度逐漸減弱。因此，從垂直速度和擾動(dòng)位溫分布可以發(fā)現(xiàn)存在明顯靜力地形波動(dòng)特點(diǎn)，波動(dòng)垂直傳播，且主要分布于山脊上方區(qū)域，位相隨高度向上游傾斜。需要說(shuō)明的是，當(dāng)參考高度選為100 m時(shí)，Eta 坐標(biāo)完全退化為Sigma 坐標(biāo)（見(jiàn)方案EXP01_），與EXP01 具有相同的模擬圖像（圖略）；伴隨參考高度逐漸減小，垂直坐標(biāo)變換為Sigma-Eta 混合坐標(biāo)，與EXP01 模擬差異逐漸顯現(xiàn)；當(dāng)參考高度為0 時(shí)，垂直坐標(biāo)為純Eta 坐標(biāo)，此時(shí)差異更加顯著（具體結(jié)果見(jiàn)4.3 節(jié)）。

圖3 EXP01 試 驗(yàn) 山 脈 背 風(fēng) 波 動(dòng) 結(jié) 構(gòu)：（a）擾 動(dòng) 位 溫（單 位：K）；（b）擾動(dòng)氣壓（單位：Pa）；（c）垂直速度（單位：m s?1）Fig.3 Structure of the mountain lee wave from experiment EXP01:(a) Potential temperature anomaly (units: K), (b) pressure anomaly(units: Pa), and (c) vertical velocity (units: m s?1)

4.3 對(duì)比試驗(yàn)結(jié)果分析

4.3.1 Eta 坐標(biāo)（參考高度50 m）

如果將參考高度下調(diào)到50 m（如圖4 所示），此時(shí)垂直速度和位溫?cái)_動(dòng)呈現(xiàn)明顯的雙中心結(jié)構(gòu)，中心強(qiáng)度明顯下降（約為EXP01 的一半）；氣壓擾動(dòng)中心強(qiáng)度亦有減小。因此，當(dāng)參考高度調(diào)整為50 m（鐘型山脈最大高度的一半）時(shí)，模擬形態(tài)保持了向上游傾斜的山脈波動(dòng)分布特征，但波動(dòng)強(qiáng)度減弱，雙中心結(jié)構(gòu)明顯。

圖4 如圖3，但為EXP02 試驗(yàn)Fig.4 Same as Fig.3 but for EXP02

4.3.2 Eta 坐標(biāo)（參考高度0）

如果取參考高度為0，則完全轉(zhuǎn)變?yōu)镋ta 坐標(biāo)。圖5 給出模擬結(jié)果，可以發(fā)現(xiàn)模擬圖像發(fā)生很大變化，基本沒(méi)有刻畫(huà)出山脈波垂直傳播、位相向上游傾斜等特征。其中，擾動(dòng)氣壓由負(fù)轉(zhuǎn)正，強(qiáng)度減弱；垂直速度和擾動(dòng)位溫強(qiáng)度非常弱，比對(duì)照試驗(yàn)減小約2～3 個(gè)數(shù)量級(jí)；另外，迎風(fēng)坡主要表現(xiàn)為弱上升氣流，背風(fēng)坡主要體現(xiàn)為弱下沉氣流。因此，Eta 垂直坐標(biāo)動(dòng)力框架在當(dāng)前模擬方案條件下，無(wú)法正確模擬靜力特性山脈波動(dòng)。

圖5 如圖3，但為EXP03 試驗(yàn)Fig.5 Same as Fig.3 but from EXP03

4.3.3 Eta 坐標(biāo)（參考高度0，垂直分辨率提高）

從以上數(shù)值試驗(yàn)分析可以發(fā)現(xiàn)，Sigma 坐標(biāo)和Eta 坐標(biāo)模式采用完全一致的數(shù)值計(jì)算方案，但對(duì)山脈波的模擬卻有很大差別：前者較好地刻畫(huà)了山脈波特征，而后者幾乎沒(méi)有模擬能力。本節(jié)將嘗試通過(guò)提高垂直分辨率來(lái)改進(jìn)Eta 坐標(biāo)WRF 新框架對(duì)山脈波的模擬性能。當(dāng)把垂直分層從80 層提高到400 層時(shí)（圖6），山脈波動(dòng)特征出現(xiàn)，重力波垂直上傳的特點(diǎn)比較明顯，位相向上游傾斜也得到體現(xiàn)。不足之處在于波動(dòng)垂直向上衰減更快，擾動(dòng)強(qiáng)度雖比EXP03 強(qiáng)，但比起對(duì)照試驗(yàn)（EXP01）仍顯著偏弱。

圖6 如圖3，但為EXP04 試驗(yàn)Fig.6 Same as Fig.3 but for EXP04

若將垂直分層繼續(xù)提高到600 層（圖略），此時(shí)山脈波模態(tài)雖仍可顯現(xiàn)，但是波動(dòng)分布較為零散，系統(tǒng)性不強(qiáng)。因此，當(dāng)使用階梯地形動(dòng)力框架模擬山脈重力波時(shí)，提高垂直分辨率可以提升模擬性能，但也不是垂直分辨率越高越好。

4.3.4 Eta 坐標(biāo)（參考高度0，模擬時(shí)長(zhǎng)加大）

我們還發(fā)現(xiàn)，在垂直分辨率較低時(shí)（如80 層），經(jīng)過(guò)長(zhǎng)時(shí)間積分調(diào)整，也能逐漸形成山脈波。從圖7 可以看出，該方案模擬出了地形波動(dòng)的垂直傳播及位相隨高度向上游傾斜的特征，具有典型靜力山脈波特征。進(jìn)一步分析發(fā)現(xiàn)，80 層垂直分層情形下出現(xiàn)典型山脈波模態(tài)的時(shí)間約為220000 s，是相同垂直分辨率Sigma 動(dòng)力框架模擬生成山脈波模態(tài)所需時(shí)間的6～7 倍。因此，對(duì)于Eta 坐標(biāo)動(dòng)力框架而言，比起相同垂直分辨率的Sigma 動(dòng)力框架需要更長(zhǎng)時(shí)間進(jìn)行波動(dòng)調(diào)整，才能在地形強(qiáng)迫下生成靜力山脈波。

圖7 如圖3，但為EXP05 結(jié)果Fig.7 Same as Fig.3 but for EXP05

5 實(shí)例試驗(yàn)

5.1 北美大陸西風(fēng)槽的模擬

5.1.1 試驗(yàn)方案

本節(jié)將利用Eta 坐標(biāo)WRF 新動(dòng)力框架對(duì)北美大陸西風(fēng)槽演變進(jìn)行模擬，并和原Sigma 坐標(biāo)框架模擬結(jié)果展開(kāi)對(duì)比。表3 給出北美大陸西風(fēng)槽演變模擬試驗(yàn)方案。當(dāng)前，模式分辨率取為30 km，時(shí)步180 s，積分24h；模式東西、南北和垂直方向分別為74、61 和28 個(gè)格點(diǎn)（對(duì)應(yīng)的地理空間范圍，水平方向：25.5°N～43.2°N，86.4°W～108°W；垂直方向：0～19 km），側(cè)邊界選為時(shí)變邊界，垂直方向剛壁；不考慮任何物理過(guò)程。

表3 北美大陸西風(fēng)槽演變模擬試驗(yàn)方案Table3 Numerical experiment scheme for simulation on the evolution of the westerly trough in the North American continent

首先分析Eta 和Sigma 兩個(gè)動(dòng)力框架生成的模式地形之異同。圖8 給出原WRF 模式地形和我們構(gòu)造的階梯地形分布。可以看出，新構(gòu)造的階梯地形與原Sigma 坐標(biāo)下追隨地形具有較好一致，如100 m 平緩地形范圍及500 m 以上較陡地形分布等。從地形高度量值看，新構(gòu)造階梯地形與原WRF 地形之差在?40 m 和160 m 之間，100 m 以上地形高度差主要分布在500 m 以上較陡地形區(qū)域。這種差別可能來(lái)自于階梯地形構(gòu)造時(shí)標(biāo)準(zhǔn)大氣狀態(tài)的選擇、Eta 垂直分層以及模式水平分辨率等因素的影響。

圖8（a）原Sigma 坐標(biāo)下和（b）Eta 坐標(biāo)下WRF 模式地形高度（單位：m）分布Fig.8 Distributions of WRF model terrain height (units: m) under (a)original Sigma coordinate core and (b) Eta coordinate core

從風(fēng)場(chǎng)的緯向—高度剖面（圖9）可以發(fā)現(xiàn)，Eta 和Sigma 模式模擬的風(fēng)場(chǎng)分布都是低層弱而高層強(qiáng)，槽區(qū)小而槽前、槽后大，且槽后更大；但Eta 坐標(biāo)新框架模擬的槽前高空風(fēng)速偏弱。圖10給出兩垂直坐標(biāo)框架下模擬24 h 的500 hPa 風(fēng)場(chǎng)分布?？梢园l(fā)現(xiàn)，選用Eta 和Sigma 垂直坐標(biāo)模擬出一致的流場(chǎng)形勢(shì)，如顯著的氣旋性流場(chǎng)分布；主要差別在Sigma 坐標(biāo)模式在（36°N，94°W）附近模擬出局地中尺度氣旋性渦旋，但Eta 坐標(biāo)模式模擬出一致的偏南氣流分布。

5.2 中國(guó)大陸槽脊系統(tǒng)的模擬

5.2.1 試驗(yàn)方案

圖9（a）Eta 坐標(biāo)新框架和（b）Sigma 坐標(biāo)原框架模擬12 h 過(guò)32°N 的水平風(fēng)速緯向—垂直剖面（單位：m s?1）Fig.9 The 12 h-simulated horizontal wind velocity across 32°N(units: m s?1): (a) Eta coordinate core; (b) original Sigma coordinate core

最后，我們來(lái)分析Eta 坐標(biāo)新框架對(duì)中國(guó)大陸低槽和副高系統(tǒng)的模擬情況。通過(guò)該方案，主要分析新框架對(duì)我國(guó)重要天氣系統(tǒng)的模擬性能。2012年7 月21 日03:00（協(xié)調(diào)世界時(shí)，下同）至22 日08:00，受貝加爾湖高空冷渦和西太平洋副熱帶高壓共同影響，華北地區(qū)出現(xiàn)了一次大范圍強(qiáng)降水天氣過(guò)程。本節(jié)主要關(guān)注相關(guān)形勢(shì)場(chǎng)的模擬情況。表4 給出模擬試驗(yàn)方案?？梢钥闯?，模擬方案與上節(jié)相同，不同之處主要體現(xiàn)在模式水平范圍，本方案涵蓋（15°N～44.7°N，109°E～144.7°E）的水平區(qū)域；另外，試驗(yàn)從2012 年7 月21 日00:00 起報(bào)，運(yùn)行12 h。

表4 中國(guó)大陸槽脊系統(tǒng)模擬試驗(yàn)方案Table4 Numerical experiment scheme for simulation on continental weather systems in China

5.2.2 試驗(yàn)結(jié)果分析

圖11 給出Eta 坐標(biāo)新框架對(duì)中國(guó)大陸低槽及副高系統(tǒng)的模擬結(jié)果，并和分析場(chǎng)進(jìn)行比較。從700 hPa 風(fēng)場(chǎng)分布可以看出，模擬和分析風(fēng)場(chǎng)分布較為一致，在副高控制區(qū)域，模擬形勢(shì)以西南風(fēng)為主，而分析則以南風(fēng)為主。

圖10（a）Eta 坐標(biāo)新框架和（b）Sigma 坐標(biāo)原框架模擬的北美大陸24 h 500 hPa 水平風(fēng)場(chǎng)（單位：m s?1）Fig.10 The 24-h simulated horizontal wind at 500 hPa (units: m s?1)in the North American continent: (a) Eta coordinate core; (b) original Sigma coordinate core

圖 11（a）Eta 坐標(biāo)新框架模擬的中國(guó)大陸6 h 700 hPa 風(fēng)場(chǎng)分布和（b）FNL 分析場(chǎng)風(fēng)場(chǎng)分布對(duì)比，單位：m s?1Fig. 11 Distributions of 6-h simulated wind at 700 hPa in the Chinese continent (units: m s?1) from (a) Eta coordinate core simulation and (b)FNL analysis

6 總結(jié)與結(jié)論

本文的主要工作是通過(guò)在WRF 模式中引入階梯地形坐標(biāo)，構(gòu)造適用于復(fù)雜地形的高精度新動(dòng)力框架；通過(guò)模式氣柱質(zhì)量變換方法，巧妙實(shí)現(xiàn)了Sigma 坐標(biāo)和Eta 坐標(biāo)下控制方程形式的完全一致，從而使Eta 坐標(biāo)WRF 新動(dòng)力框架計(jì)算更加便利、高效。

Eta 坐標(biāo)下WRF 動(dòng)力框架搭建完成后，設(shè)計(jì)理想和實(shí)例試驗(yàn)對(duì)新框架進(jìn)行檢驗(yàn)。主要結(jié)論如下：

（1）從動(dòng)力框架形式來(lái)看，Eta 垂直坐標(biāo)和Sigma 垂直坐標(biāo)下WRF 的差別主要體現(xiàn)在模式頂與參考面之間氣柱的質(zhì)量定義，前者為P2=(πs?πt)/ηs， 后者為 μ=πs?πt。通過(guò)該質(zhì)量變換設(shè)計(jì)，使得兩個(gè)動(dòng)力框架形式完全一致，從而使程序?qū)崿F(xiàn)變得簡(jiǎn)單。

（2）引入AREM 參考大氣狀態(tài)用于構(gòu)造WRF新動(dòng)力框架階梯地形，雖然該大氣狀態(tài)仍有改進(jìn)空間，但其主要用于計(jì)算地表Eta（ ηs）的初估值，對(duì)階梯地形設(shè)計(jì)和模式積分計(jì)算影響不大。此外，通過(guò)設(shè)計(jì)階梯地形2D 和3D 網(wǎng)格表征量使模式循環(huán)迭代與高階差分計(jì)算過(guò)程更加簡(jiǎn)潔、高效。

（3）理想和實(shí)例試驗(yàn)表明，Eta 坐標(biāo)WRF 動(dòng)力框架運(yùn)行穩(wěn)定、計(jì)算合理。階梯地形坐標(biāo)動(dòng)力框架可對(duì)靜力地形波進(jìn)行有效刻畫(huà)，但前提是提高垂直分辨率或延長(zhǎng)模擬時(shí)間。再考慮階梯地形先判斷再計(jì)算的額外時(shí)耗，相比Sigma 坐標(biāo)動(dòng)力框架，新框架計(jì)算花費(fèi)更大。

本文主要圍繞Eta 垂直坐標(biāo)WRF 動(dòng)力框架設(shè)計(jì)、實(shí)現(xiàn)及其驗(yàn)證展開(kāi)，數(shù)值試驗(yàn)尤其是實(shí)例試驗(yàn)分析仍顯初步，Eta 和Sigma 垂直坐標(biāo)的差異及其原因有待今后深入分析討論。另外，模式現(xiàn)有初始化僅由插值實(shí)現(xiàn)，物理過(guò)程參數(shù)化也未接入框架，也不利于實(shí)例試驗(yàn)開(kāi)展。