中国气象学会主办。
文章信息
- 李 超, 陈德辉, 李兴良. 2012.
- LI Chao, CHEN Dehui, LI Xingliang. 2012.
- 大气数值预报模式中高度地形追随坐标的设计研究:理论分析与理想试验
- A design of height-based terrain-following coordinates in the atmospheric numerical model:Theoretical analysis and idealized tests
- 气象学报, 70(6): 1247-1259
- Acta Meteorologica Sinica, 70(6): 1247-1259.
- http://dx.doi.org/10.11676/qxxb2012.105
-
文章历史
- 收稿日期:2011-04-19
- 改回日期:2012-03-07
2. 中国科学院大学,北京,100049;
3. 中国气象局数值预报中心,北京,100081
2. University of Chinese Academy of Science, Beijing 100049, China;
3. Center of Numerical Weather Prediction, China Meteorological Administration, Beijing 100081, China
随着计算机性能的提高,高分辨率的非静力数值预报模式得到了迅速发展。例如,中国气象局的新一代全球区域同化预报多尺度统一模式GRAPES、美国中尺度模式WRF、德国LM模式等。在非静力模式的设计中,选取适当的垂直坐标描述非静力特性非常重要;同时,模式水平分辨率提高使得可分辨的地形尺度减小、地形坡度增大,选取的垂直坐标既要能描述地形影响,又要能够减小地形作用带来的计算误差显得很关键。
自数值天气预报系统产生以来,垂直坐标的发展经历了由“非归一化”坐标向“归一化”坐标的转变,这种转变主要是为了应对下边界处理中遇到的困难以及适应非静力模式的发展需求。“非归一化”垂直坐标主要有几何高度坐标(Richardson,1965)、气压坐标(Eliassen,1949)以及适用于锋面研究的位温坐标。但是,这3种坐标在地形起伏的地区很难给出准确的下边界条件。“归一化”垂直坐标主要有气压地形追随坐标和高度地形追随坐标。Phillips(1957)首先提出了σ坐标系(即气压地形追随坐标),这种对气压归一化的坐标便于计算机在整个计算域内进行运算,地形层上的垂直风经坐标转换后消失,简化了下边界条件,同时可变化间距的计算层为模式动力与各种边界层的相关参数化方案的耦合提供了便利。Kasahara(1974)将σ坐标的特点引伸到几何高度坐标和位温坐标中,使用通用的坐标形式给出了静力大气下水平运动方程、静力平衡方程、质量守恒方程和热力学方程等一般表达式。为了描述大气的非静力特征,Gal-Chen等(1975)提出了高度地形追随坐标。Schar等(2002)提出了平缓坐标(SLEVE,smooth level vertical coordinate),考虑不同尺度地形的影响,在高层使坐标面更加平缓。李兴良等(2005)也分析比较了非静力模式中的气压地形追随坐标和高度地形追随坐标,发现高度地形追随坐标具备一定的预报优势。
尽管垂直坐标发展为“归一化”坐标解决了复杂的下边界问题,但也带来了新的问题。气压梯度力的计算由一项变成了两个大项的小差,容易产生计算误差,尤其在地形复杂陡峭处这种误差更加明显,因而很大程度上影响数值预报的结果。为了解决这个问题,主要是从以下几个方向研究的:第一是不改变垂直坐标形式对气压梯度力项的计算方案进行改进,例如:回插等压面法、递推算法(杨晓娟等,2003)、扣除法(钱永甫等,1994;周天军等,1996)、特殊差分格式法、差微差算法(钱永甫等,1991)、气压回插法(Hu et al,2007)。第二是特别处理地形,典型例子就是η阶梯坐标(Mesinger,1984)的提出,η坐标将模式地形处理为“台阶地形”,简化了下边界条件。第三是为改变坐标面随高度的衰减状况提出新的地形追随坐标,将垂直坐标的选择转变为衰减系数的选择(Schar et al,2002)。第三种方法的优点是通过对衰减系数的调节,使垂直坐标在底层受到地形作用的影响,而在高层坐标面趋于平直,更接近大气的真实状况,这也是本文研究的出发点。
Gal-Chen&Somerville(Gal.C.S)坐标、单尺度平缓坐标(SLEVE1)、双尺度平缓坐标(SLEVE2)以及三角函数平缓坐标(COS)都是比较典型的高度地形追随坐标,其性质因衰减系数而不同。本研究对这几种坐标进行了分析,比较了在有地形和无地形情况下不同坐标对气压梯度力计算误差以及质量平流输送的影响,为模式垂直坐标的选择提供参考。以下如无特别说明,所提到的坐标均为高度地形追随坐标。 2 高度地形追随坐标理论分析
Gal-Chen等(1975)指出,地形追随坐标的设计必须满足两个条件:(1)坐标转换前后两者的计算域必须一致,即:在z=ZS处,=0;在z=ZT处,=ZT;(2)坐标转换必须是单调可逆的,即坐标转换的Jb-1项必须满足。
Jb-1是坐标转换过程中的重要参数,在下面的理论分析中将重点讨论这个参数。 2.1 几种高度地形追随坐标
任何一种高度地形追随坐标可以写成统一形式
(1)Gal.C.S坐标的定义式为 =,将其对应式(1)转化为z=+(1-)ZS,可以看出衰减系数为
由式(2)可知,衰减系数b随模式高度的衰减率为常数-。
(2)SLEVE平缓坐标将地形分离为多种尺度,对不同尺度的地形选择不同的衰减系数。SLEVE1坐标地形未划分尺度,其衰减系数为
(3)SLEVE双尺度平缓坐标。事实上,模式地形高度的差异很大,高可达数千米(如珠穆朗玛峰),低则只有几十米(如丘陵),几乎可以用连续尺度谱ZSi(i=1,2,3…,L)来表示,对于不同的模式地形高度“尺度谱段”,式(3)中的h*也应选择不同的地形高度特征尺度h*i,式(3)可以改写为更一般的形式
(4)COS坐标。SLEVE坐标是以“双曲函数”为衰减函数基的高度地形追随坐标,这里选择“三角函数”为衰减函数基构造高度地形追随坐标(Klemp,2011),即
衰减系数b将随模式层高度升高而减小,其使坐标面随增大由曲面向平面衰减。当b=0时,模式面变成了水平面,取不同的b有不同的衰减影响。这样对垂直坐标的改进和选择问题就转化成了对衰减系数的改进和选择问题。在统一的坐标形式下,只对衰减系数b进行改变,程序实现更加方便简洁。衰减系数b成为坐标转换过程中的重要参数,下面分析中将着重讨论这个参数。图 1是衰减系数b各不相同的4种坐标的坐标面分布情况。
2.2 几种坐标的性质分析可以较容易证明在4种坐标中z=ZS处,=0,在z=ZT处,=ZT,Gal.C.S和SLEVE1坐标满足坐标转换的单调性条件为Jb-1>0,对于SLEVE2和COS坐标还需要一些约束条件才能满足Jb-1>0。
SLEVE2坐标的坐标转换项为
COS坐标的转换项为
由式(9)可知,Jb-1受到ZS、Zc、n的共同约束。例如,模式层顶高度ZT为25 km,临界衰减高度Zc为10 km时,n若为1,则Jb-1在临界衰减高度处不连续,n取值要大于1。如图 2所示,n≤12时,Jb-1>0;n≥13(虚线为n=13)时,Jb-1<0,则n取值不可大于12。倘若模式层顶高度、地形高度、临界衰减高度等发生变化,满足Jb-1>0的n的最大值也会发生变化。下面二维质
量平流试验中n取值6时,设置临界衰减高度分别为10和15 km均可满足Jb-1>0,其对应的最大地形高度分别约为4.1和6.1 km。
2.3 不同坐标地形高度衰减率的比较由式(1)可知,不同坐标的差别在于对地形作用的衰减效果的不同,即bZS的差别。下面选择地形坡度最大的点比较4种坐标的地形高度衰减率∂[bZS(x,y)]/∂z。首先从地形高度衰减率的角度对SLEVE2坐标进行理论分析,其次比较COS坐标中n取不同值时地形高度衰减率的分布,最后选择合适的n值比较4种坐标的地形高度衰减率。
地形函数满足
SLEVE2坐标选取的大尺度和小尺度地形分别为
图 3a为从SLEVE2坐标中分离出的大、小尺度地形分布。根据实践和坐标转换的单调性要求分别取大尺度和小尺度地形的衰减高度h*1和h2*为15和2.5 km,对应式(4)中的衰减系数,这样大尺度地形和小尺度地形的衰减状况将有所不同。图 3b为实际地形、大尺度和小尺度地形的地形高度衰减率分布。可以看出,在15 km以下小尺度
地形的衰减远快于大尺度地形,在低层,实际地形衰减主要体现在小尺度地形的衰减上。
COS坐标中涉及到参数n的选择,从COS坐标中参数n分别取2、6、8、12时,∂(bZS)/∂z的比较(图 4)可以看出,COS坐标在临界衰减高度以下将按地形作用进行衰减,且n取不同值时最大衰减率出现的高度不同,n越大这个高度越低,地形衰减相对就越快。
下面分析Gal.C.S、SLEVE1、SLEVE2和COS坐标(Zc=10 km,n=6)的地形高度衰减率。图 5a为4种坐标下各层坐标面的坡度分布。Gal.C.S坐标各层坐标面坡度线性递减;SLEVE2坐标的坐标面坡度比Gal.C.S坐标和SLEVE1坐标减小得快,15 km时SLEVE2坐标的坐标面坡度已经减小到很小;COS坐标在大约4 km以上坐标面坡度小于SLEVE2坐标,4 km以下则大于SLEVE2坐标。图 5b为4种坐标下地形高度衰减率∂(bZS)/∂z 的分布。Gal.C.S坐标地形高度衰减率为常值;在同一高度上SLEVE2坐标和COS坐标的地形高度衰减率比其他两种坐标都要大,对地形作用的衰减相对快一些;COS坐标在大约4 km出现地形高度衰减率的最大值,在这一高度以上比SLEVE2坐标衰减快,但这一高度以下则相对较慢。可见,地形高度衰减率的分布与各层坐标面的坡度有很好的对应关系,地形高度衰减率直接影响坐标面的平滑程度。
从理论分析的结果可以看出,SLEVE2坐标将地形区分尺度考虑,把尺度较小的地形迅速衰减掉。与地形滤波相比,在保持地形最大高度不变的前提下,降低了可识别的地形坡度的量级,比只考虑单一尺度地形的SLEVE1坐标有优势。COS坐标给出了连续的衰减系数,在高层随高度增加对地形的衰减比SLEVE2坐标更快。3 不同坐标气压梯度力计算误差分析
“归一化”坐标带来的主要问题是气压梯度力的计算误差变大,特别是在地形起伏较大的地区。下面比较4种坐标下气压梯度力的计算对地形坡度的敏感性。
试验比较在有地形的静止大气中,4种坐标气压梯度力的计算误差的差异。基于理论分析的结果,COS坐标中参数设置为Zc=10 km,n=6,地形分布满足式(10),参考大气取为
与图 6气压梯度力误差直观图相比,下面将其误差量化,选择气压梯度力误差最大的点,依据相对误差计算公式
比较SLEVE2和COS坐标在第20层(9.5 km)以下的气压梯度力误差分布(图 7)可知,COS坐标在各个层上气压梯度力误差的分布和地形坡度有很好的对应,地形坡度大(小)之处气压梯度力计算误差也大(小),但SLEVE2坐标中却不存在这种对应关系,这是因为SLEVE2坐标对不同尺度谱段的地形高度的衰减不同。在实践中可以根据需要对有可能产生大的气压梯度力计算误差的地形影响进行快速衰减。在第10层(4.5 km)以下,COS比SLEVE2坐标的误差要大一些,而在第10层以上比SLEVE2坐标小。
综上所述,SLEVE2、COS坐标(Zc=10 km,n=6)与Gal.C.S、SLEVE1坐标相比,可以更好地减小坐标面上的气压梯度力计算误差。SLEVE2坐标可以解决地形精细化引发的问题,COS坐标可以在一定限度内让地形影响迅速衰减,两种坐标各具优势。4 不同坐标下的二维质量平流试验
由上述理论分析可知,大气数值预报模式中不同坐标对地形作用的衰减性质不同。二维质量平流试验可较好地检验大气数值预报模式中采用上述不同高度地形追随坐标时地形对物质输送的影响。该试验采用密度通量方程,坐标转换后其通用表达式为
为了计算水平风速u、v,需要在方程中引入流函数来估算速度场,设流函数为φ(z)=-∫0zu(z)dz,由坐标转换后,气团的平流速度改写为,整理后密度通量方程转化为
试验设计如图 8所示,有一气团初始时刻(t=0 s)密 度ρ分布满足
半径宽Rx=25 km,高Rz=3 km,气团位于(x0,z0)=(-50 km,9 km)。无误差的情况下:平流运动5000 s后,气团位于山头之上,即(x2,z2)=(0 km,9 km);平流运动10000 s后,气团位于山头的另一侧,与初始时刻的气团位置相对称,即(x3,z3)=(50 km,9 km)。图 8为平流试验的密度解析值分布,给定气团的平流速度为
质量通量方程的差分计算采用Arakawa-c网格(Arakawa et al,1977)和有限体积的一阶迎风方案。令式(16)中,对通量项采取分维算法
由上述分析可知,在有地形情况下,不同坐标的坐标面起伏程度(亦称坡度)不同(图 1),因而对质量输送的影响也不相同,可以通过比较输送过程中某时刻的密度场及其变化状况来说明不同坐标对质量输送的影响。根据Williamson等(1992)提出的计算全球积分误差的方法,本研究定义衡量误差公式为
下面分别比较在有地形和无地形情况下,t=0、5000、10000 s时4种坐标的密度分布,t=10000 s时的质量平流误差分布,以及积分各个时刻质量平流标准化误差演变。
图 9是一阶迎风方案下的4种坐标质量平流试验结果,其中图 9a6、b6是无地形数值模拟结果,可作为该数值计算方案的参考解。图 9a1—a6为不同坐标下对应于图 8的3个时刻的密度场模拟结果。与图 8的解析值比较可以看出,一阶迎风方案下无地形时质量气团随时间积分出现了水平拉伸,这主要是数值计算方案的耗散性决定的。COS坐标(Zc=10 km)与无地形结果最接近,这是在临界衰减高度Zc=10 km以下该坐标对地形平滑的结果;SLEVE2与COS坐标(Zc=15 km)在t=5000 s时气团形变相对于Gal.C.S与SLEVE1坐标都较小,只在气团的下部有一些地形影响特征;Gal.C.S与SLEVE1坐标在经过地形上方时气团形变较严重,并带有明显的地形影响特征,t=10000 s时气团形状与无地形时相差甚远。图 9b1—b6为t=10000 s的质量平流误差(即该时刻密度数值解与解析解之差)分布。4种坐标与无地形时在该时刻有不同的质量平流误差,且误差从Gal.C.S、SLEVE1、SLEVE2到COS坐标逐渐减小;无地形情况下绝对误差最大值为0.22,Gal.C.S、SLEVE1、SLEVE2坐标绝对误差最大值分别为0.65、0.61、0.38,而COS坐标(Zc=10,15 km)与无地形最接近,误差最大值为0.32;特别是COS(Zc=15 km)在较高层比SLEVE2坐标绝对误差小。
图 10是4种坐标和无地形情况下标准化误差的时间演变。图 10a为l2误差演变,可以看出,4种坐标下平流输送的l2标准化误差均随时间增大,在2000—7000 s误差迅速增大,这是因为2000—7000 s气团恰好通过地形上方,可以看出,地形对平流误差的影响较大;而7000 s以后误差变化不大。相比较而言,SLEVE2和COS坐标的误差是4种坐标中较小的,且COS坐标Zc=10 km比Zc=15 km时误差小,Zc=10 km时误差最接近无地形情况。图 10b为误差演变l∞情况,其误差分布与图 10a大致类似,它反映了最大误差(极值)随时间的演变。特别是COS坐标Zc=10 km时误差与无地形时基本重合。不论是l2还是l∞,10000 s时Gal.C.S坐标的标准化误差比无地形大50%左右,而SLEVE2和COS坐标相对更接近无地形情况。
由试验结果和误差分析可以看出,平流误差与坐标面的起伏有关。SLEVE2、COS与Gal.C.S、SLEVE1坐标相比越往高层坐标面越平缓,因而可以更好地减小平流误差;并且,COS坐标临界衰减高度Zc越小越接近无地形时的误差。5 结论和讨论
选取了4种典型的高度地形追随坐标,将其化为统一的坐标形式,分析其地形衰减系数b、坐标转换的雅可比项Jb以及坐标面上气压梯度力计算误差,然后分别进行了二维质量平流试验。通过理论分析和理想试验,可以得到以下几个初步结论:
(1)高度地形追随坐标可以转化为统一形式z=+bZS,将对坐标的研究简化为对地形衰减系数b的研究。
(2)SLEVE2坐标要满足坐标转换的Jb-1>0条件,需要根据实践和理论合理地组合地形高度的各个尺度谱段ZSi及与其对应的特征尺度h*i(i=1,2)。
(3)COS坐标的Jb-1受到地形高度ZS、衰减临界高度Zc和参数n的约束。要满足坐标转换的Jb-1>0条件:ZS和Zc确定时,n取值越大坐标面衰减效果越好,n取值有上界,且上界范围内n也不宜取值过大;Zc和n确定时,地形高度也有上界,若取临界衰减高度Zc为10和15 km,地形最大高度分别约为4.1和6.1 km。
(4)理论分析Gal.C.S、SLEVE1、SLEVE2和COS坐标可以发现,SLEVE2和COS坐标的衰减效果要明显好于其他2种坐标;SLEVE2和COS坐标相比,高层COS坐标衰减快,低层SLEVE2坐标衰减快。
(5)气压梯度力误差分析表明,在同一模式层上,SLEVE2和COS坐标比Gal.C.S和SLEVE1坐标的气压梯度力计算误差小;在某一高度以上COS坐标气压梯度力计算误差小于SLEVE2坐标,这一高度以下情况相反。
(6)二维理想平流试验结果表明,SLEVE2和COS坐标减小平流误差效果优于Gal.C.S和SLEVE1坐标;COS坐标中临界衰减高度Zc的选择很重要,Zc越小越接近无地形的结果。
以上讨论的结果是基于理论分析和理想数值试验,各种高度地形追随坐标在真实大气数值预报模式如GRAPES模式(张人禾等,2008;陈德辉等,2008;杨学胜等,2008)中的表现还需要实际检验;同时,可以考虑结合各种坐标的优点设计新的高度地形追随坐标。
陈德辉, 薛纪善, 杨学胜等. 2008. GRAPES新一代全球/区域多尺度统一数值预报模式总体设计研究. 科学通报, 53(20): 2396-2407 |
李兴良, 陈德辉. 2005. 非静力中尺度高分辨率模式模拟中的垂直坐标影响研究. 气象学报, 63(2): 161-171 |
钱永甫, 王云峰. 1991. 数值模式中气压梯度力的算法试验. 气象学报, 49(3): 538-547 |
钱永甫, 周天军. 1995. 有地形模式中气压梯度力误差扣除法. 高原气象, 14(1): 1-9 |
钱永甫, 周天军. 1994. 陡峭地形区气压梯度力的误差扣除法. 热带气象学报, 10(4): 358-368 |
杨晓娟, 钱永甫. 2003. P-σ坐标系数值模式中气压梯度力的递推算法. 大气科学, 27(2): 171-181 |
杨学胜, 胡江林, 陈德辉等. 2008. 全球有限区数值预报模式动力框架的试验验证. 科学通报, 53(20): 2418-2424 |
张人禾, 沈学顺. 2008. 中国国家级新一代业务数值预报系统GRAPES的发展. 科学通报, 53(20): 2393-2395 |
周天军, 钱永甫. 1996. 地形区两种典型气压梯度力计算方法的比较. 应用气象学报, 7(3): 377-380 |
Arakawa A, Lamb V R. 1977. Computational design of the basic dynamical processes of the UCLA general circulation model//Methods in Computational Physics, General Circulation Models of the Atmosphere. New York:Academic Press, 17 |
Eliassen A. 1949. The quasi-static equations of motion with pressure as independent variable.Geofys Publikasjoner, 17(3):44p |
Gal-Chen T J, Somerville R C. 1975. On the use of a coordinate transformation for the solution of the Navier-Stokes equations. J Comput Phys, 17(2): 209-228 |
Hu J L, Wang P X. 2007. The errors of pressure gradient force in high-resolution Meso-scale model with terrain-following coordinate and its revised scheme. J Atmos Sci, 31(1): 110-118 |
Kasahara A. 1974. Various vertical coordinate systems used for numerical weather prediction. Mon Wea Rev, 102(7): 509-522 |
Klemp J B. 2011. A terrain-following coordinate with smoothed coordinate surfaces. Mon Wea Rev, 139(7): 2163-2169 |
Mesinger F. 1984. A blocking technique for representation of mountains in atmospheric models. Riv Meteor Aeronautica, 44(2): 195-20 |
Phillips N A. 1957. A coordinate system having some special advantages for numerical forecasting. J Meteor, 14(2): 184-185 |
Richardson L F. 1965. Weather Prediction by Numerical Process. New York: Cambridge University Press, 236pp |
Schar C, Leuenberger D, Fuhrer O, et al. 2002. A new terrain-following vertical coordinate formulation for atmospheric prediction models. Mon Wea Rev, 130(10): 2459-2480 |
Williamson D L, Drake J B, Hack J J, et al. 1992. A standard test set for numerical approximations to the shallow water equations in spherical geometry. J Comput Phys, 102(1): 211-224 |