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

1 地理信息工程國家重點實驗室，西安710054

2 中國氣象局，北京 100081

3 中國科學院大氣物理研究所大氣科學和地球流體力學數(shù)值模擬國家重點實驗室(LASG)，北京 100029

1 引言

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

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

1.1 垂直坐標設(shè)計

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

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

1.2 地形追隨坐標的問題及其改進

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

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

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

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

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

1.3 山脈波的求解

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

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

2 Eta 坐標下WRF 動力框架

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

（a）第一運動方程：

（b）第二運動方程：

（c）第三運動方程：

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

（e）熱力學方程：

（f）擾動位勢方程：

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

（i）靜力近似方程：

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

其中，

式中， π表示靜力氣壓， πt為模式頂靜力氣壓， πs為地面靜力氣壓， ηs表 示地形處Eta 值； πrf為參考靜力氣壓；zs為 地形高度，zb為 平緩地形高度， α表示大氣比容，p0表示參考氣壓，一般取1000 hPa。

其他變量使用氣象常用符號表示，不再贅述。

3 Eta 垂直坐標設(shè)計

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

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

3.2 階梯地形表征

3.2.1 階梯地形2D 表征量

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

圖1（a）Eta-WRF 模式的階梯地形高度（單位：m）及階梯坐標模式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

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

3.2.2 3D 表征量

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

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

4 山脈波試驗

4.1 試驗方案

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

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

表2 Eta 坐標對比試驗方案Table2 Configuration for comparative experiments with Eta coordinates

4.2 對照試驗結(jié)果分析

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

圖3 EXP01 試 驗 山 脈 背 風 波 動 結(jié) 構(gòu)：（a）擾 動 位 溫（單 位：K）；（b）擾動氣壓（單位：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 對比試驗結(jié)果分析

4.3.1 Eta 坐標（參考高度50 m）

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

圖4 如圖3，但為EXP02 試驗Fig.4 Same as Fig.3 but for EXP02

4.3.2 Eta 坐標（參考高度0）

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

圖5 如圖3，但為EXP03 試驗Fig.5 Same as Fig.3 but from EXP03

4.3.3 Eta 坐標（參考高度0，垂直分辨率提高）

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

圖6 如圖3，但為EXP04 試驗Fig.6 Same as Fig.3 but for EXP04

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

4.3.4 Eta 坐標（參考高度0，模擬時長加大）

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

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

5 實例試驗

5.1 北美大陸西風槽的模擬

5.1.1 試驗方案

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

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

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

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

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

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

5.2.1 試驗方案

圖9（a）Eta 坐標新框架和（b）Sigma 坐標原框架模擬12 h 過32°N 的水平風速緯向—垂直剖面（單位：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

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

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

5.2.2 試驗結(jié)果分析

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

圖10（a）Eta 坐標新框架和（b）Sigma 坐標原框架模擬的北美大陸24 h 500 hPa 水平風場（單位：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 坐標新框架模擬的中國大陸6 h 700 hPa 風場分布和（b）FNL 分析場風場分布對比，單位：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é)論

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

Eta 坐標下WRF 動力框架搭建完成后，設(shè)計理想和實例試驗對新框架進行檢驗。主要結(jié)論如下：

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

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

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

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