1 引 言

随着计算机性能的提高，高分辨率的非静力数值预报模式得到了迅速发展。例如，中国气象局的新一代全球区域同化预报多尺度统一模式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=Z_S处，=0；在z=Z_T处，=Z_T；(2)坐标转换必须是单调可逆的，即坐标转换的J_b^-1项必须满足。

J_b^-1是坐标转换过程中的重要参数，在下面的理论分析中将重点讨论这个参数。 2.1 几种高度地形追随坐标

任何一种高度地形追随坐标可以写成统一形式

式中，为高度地形追随坐标的高度(即模式层高度)，z为自然几何高度，Z_S为模式地形高度。由式(1)可以看出，b即是高度地形追随坐标的衰减系数。下面给出几种典型坐标的b的表达式。

(1)Gal.C.S坐标的定义式为 =，将其对应式(1)转化为z=+(1-)Z_S，可以看出衰减系数为

式中，Z_T为模式层顶高度。

由式(2)可知，衰减系数b随模式高度的衰减率为常数-。

(2)SLEVE平缓坐标将地形分离为多种尺度，对不同尺度的地形选择不同的衰减系数。SLEVE1坐标地形未划分尺度，其衰减系数为

式中，h^*为模式地形高度的特征尺度，模式坐标面均以这一特征尺度的地形高度进行衰减。这里的h^*相当于Schar等(2002)中的s。

(3)SLEVE双尺度平缓坐标。事实上，模式地形高度的差异很大，高可达数千米(如珠穆朗玛峰)，低则只有几十米(如丘陵)，几乎可以用连续尺度谱ZS_i(i=1，2，3…，L)来表示，对于不同的模式地形高度“尺度谱段”，式(3)中的h^*也应选择不同的地形高度特征尺度h^*_i，式(3)可以改写为更一般的形式

其对应的坐标转换关系式(1)则改写为该坐标称为SLEVE多尺度平缓坐标。若取L=2，这时的坐标称为SLEVE双尺度平缓坐标(SLEVE2坐标)，该坐标将地形划分为大尺度地形Z_h^*1和小尺度地形Z_h2^*，其对应的衰减系数分别为

(4)COS坐标。SLEVE坐标是以“双曲函数”为衰减函数基的高度地形追随坐标，这里选择“三角函数”为衰减函数基构造高度地形追随坐标(Klemp，2011)，即

式中，Z_c是临界衰减高度，n是正整数。

衰减系数b将随模式层高度升高而减小，其使坐标面随增大由曲面向平面衰减。当b=0时，模式面变成了水平面，取不同的b有不同的衰减影响。这样对垂直坐标的改进和选择问题就转化成了对衰减系数的改进和选择问题。在统一的坐标形式下，只对衰减系数b进行改变，程序实现更加方便简洁。衰减系数b成为坐标转换过程中的重要参数，下面分析中将着重讨论这个参数。图 1是衰减系数b各不相同的4种坐标的坐标面分布情况。

2.2 几种坐标的性质分析

可以较容易证明在4种坐标中z=Z_S处，=0，在z=Z_T处，=Z_T，Gal.C.S和SLEVE1坐标满足坐标转换的单调性条件为J_b^-1>0，对于SLEVE2和COS坐标还需要一些约束条件才能满足J_b^-1>0。

SLEVE2坐标的坐标转换项为

这是由分离出的两种尺度的地形高度Z_h^*1和Z_h2^*及它们对应的衰减高度h^*₁和h₂^*共同约束的。不等式是上式成立的必要但不充分条件，其中，Z_h^*1，max和Z_h2^*，max分别为大尺度地形和小尺度地形的最大高度。可见只需要找到的合理组合即可满足J_b^-1>0。实践上，h^*₁一般取值为10 km左右，然后再基于此给出使小尺度地形迅速衰减的h₂^*的值。

COS坐标的转换项为

由式(9)可知，J_b^-1受到Z_S、Z_c、n的共同约束。例如，模式层顶高度Z_T为25 km，临界衰减高度Z_c为10 km时，n若为1，则J_b^-1在临界衰减高度处不连续，n取值要大于1。如图 2所示，n≤12时，J_b^-1＞0；n≥13(虚线为n=13)时，J_b^-1＜0，则n取值不可大于12。倘若模式层顶高度、地形高度、临界衰减高度等发生变化，满足J_b^-1＞0的n的最大值也会发生变化。下面二维质

量平流试验中n取值6时，设置临界衰减高度分别为10和15 km均可满足J_b^-1＞0，其对应的最大地形高度分别约为4.1和6.1 km。

2.3 不同坐标地形高度衰减率的比较

由式(1)可知，不同坐标的差别在于对地形作用的衰减效果的不同，即bZ_S的差别。下面选择地形坡度最大的点比较4种坐标的地形高度衰减率∂[bZ_S(x，y)]/∂z。首先从地形高度衰减率的角度对SLEVE2坐标进行理论分析，其次比较COS坐标中n取不同值时地形高度衰减率的分布，最后选择合适的n值比较4种坐标的地形高度衰减率。

地形函数满足

式中，3 km，a=2.5 km，λ_S=8 km，最大地形坡度约为1.4698。

SLEVE2坐标选取的大尺度和小尺度地形分别为

图 3a为从SLEVE2坐标中分离出的大、小尺度地形分布。根据实践和坐标转换的单调性要求分别取大尺度和小尺度地形的衰减高度h^*₁和h₂^*为15和2.5 km，对应式(4)中的衰减系数，这样大尺度地形和小尺度地形的衰减状况将有所不同。图 3b为实际地形、大尺度和小尺度地形的地形高度衰减率分布。可以看出，在15 km以下小尺度

地形的衰减远快于大尺度地形，在低层，实际地形衰减主要体现在小尺度地形的衰减上。

COS坐标中涉及到参数n的选择，从COS坐标中参数n分别取2、6、8、12时，∂(bZ_S)/∂z的比较(图 4)可以看出，COS坐标在临界衰减高度以下将按地形作用进行衰减，且n取不同值时最大衰减率出现的高度不同，n越大这个高度越低，地形衰减相对就越快。

下面分析Gal.C.S、SLEVE1、SLEVE2和COS坐标(Z_c=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种坐标下地形高度衰减率∂(bZ_S)/∂z 的分布。Gal.C.S坐标地形高度衰减率为常值；在同一高度上SLEVE2坐标和COS坐标的地形高度衰减率比其他两种坐标都要大，对地形作用的衰减相对快一些；COS坐标在大约4 km出现地形高度衰减率的最大值，在这一高度以上比SLEVE2坐标衰减快，但这一高度以下则相对较慢。可见，地形高度衰减率的分布与各层坐标面的坡度有很好的对应关系，地形高度衰减率直接影响坐标面的平滑程度。

从理论分析的结果可以看出，SLEVE2坐标将地形区分尺度考虑，把尺度较小的地形迅速衰减掉。与地形滤波相比，在保持地形最大高度不变的前提下，降低了可识别的地形坡度的量级，比只考虑单一尺度地形的SLEVE1坐标有优势。COS坐标给出了连续的衰减系数，在高层随高度增加对地形的衰减比SLEVE2坐标更快。3 不同坐标气压梯度力计算误差分析

“归一化”坐标带来的主要问题是气压梯度力的计算误差变大，特别是在地形起伏较大的地区。下面比较4种坐标下气压梯度力的计算对地形坡度的敏感性。

试验比较在有地形的静止大气中，4种坐标气压梯度力的计算误差的差异。基于理论分析的结果，COS坐标中参数设置为Z_c=10 km，n=6，地形分布满足式(10)，参考大气取为

式中，Π和θ为z的函数。可以证明，在参考大气下气压梯度力为式(13)的离散采用气压梯度力计算的经典法，即对 温度计算格点的值，其余各项采用中央差分(钱永甫等，1995)。试验模式层顶高度为25 km，垂直均匀分50层，第1层为地表面；水平方向分辨率为1 km。从第2、10、20、30、40模式层(即0.5、4.5、9.5、14.5、19.5 km)上气压梯度力计算误差(图 6)可以看出，在相同的地形和计算差分格式下，同一坐标在不同模式层上引起的气压梯度力误差不同，不同坐标在相同的模式层上引起的气压梯度力误差不同。在纵向上，4种坐标的气压梯度力误差均随高度减小，这与图 5a坐标面坡度的变化相符；4种坐标气压梯度力计算误差的减小速度是不一样的，SLEVE2和COS坐标比Gal.C.S和SLEVE1坐标减小快得多。在横向上，第2层(0.5 km)SLEVE2比其余3种坐标的误差小，第10层(4.5 km)以上SLEVE2和COS坐标比其他2种坐标误差小很多；在第20层以上的某一高度上COS坐标的误差减小为0。

与图 6气压梯度力误差直观图相比，下面将其误差量化，选择气压梯度力误差最大的点，依据相对误差计算公式

式中，EGal是Gal.C.S坐标的误差，ESLEVE1/SLEVE2/COS是SLEVE1、SLEVE2或COS坐标的误差。比较该点处的第2、10、20、30、40模式层上SLEVE1、SLEVE2、COS坐标相对于Gal.C.S坐标的误差减小量。如表 1所示，在4.5 km(10层)以上，SLEVE2、COS相对Gal.C.S坐标误差减小率迅速增大到75%以上，远大于SLEVE1坐标的31%；在0.5 km(2层)处与Gal.C.S相比，COS、SLEVE1坐标误差减小的效果不明显，仅为2%和4%，但 14.5 km(30层)以上，COS坐标气压梯度力误差减小率为100%，SLEVE2坐标也达到了99%。

比较SLEVE2和COS坐标在第20层(9.5 km)以下的气压梯度力误差分布(图 7)可知，COS坐标在各个层上气压梯度力误差的分布和地形坡度有很好的对应，地形坡度大(小)之处气压梯度力计算误差也大(小)，但SLEVE2坐标中却不存在这种对应关系，这是因为SLEVE2坐标对不同尺度谱段的地形高度的衰减不同。在实践中可以根据需要对有可能产生大的气压梯度力计算误差的地形影响进行快速衰减。在第10层(4.5 km)以下，COS比SLEVE2坐标的误差要大一些，而在第10层以上比SLEVE2坐标小。

综上所述，SLEVE2、COS坐标(Z_c=10 km，n=6)与Gal.C.S、SLEVE1坐标相比，可以更好地减小坐标面上的气压梯度力计算误差。SLEVE2坐标可以解决地形精细化引发的问题，COS坐标可以在一定限度内让地形影响迅速衰减，两种坐标各具优势。4 不同坐标下的二维质量平流试验

由上述理论分析可知，大气数值预报模式中不同坐标对地形作用的衰减性质不同。二维质量平流试验可较好地检验大气数值预报模式中采用上述不同高度地形追随坐标时地形对物质输送的影响。该试验采用密度通量方程，坐标转换后其通用表达式为

为了计算水平风速u、v，需要在方程中引入流函数来估算速度场，设流函数为φ(z)=-∫₀^zu(z)dz，由坐标转换后，气团的平流速度改写为，整理后密度通量方程转化为

引入流函数后，该式可以较好地描述坐标面上的无辐散运动，可以看到J_b^-1在通量项中消失，这样避免了计算。 4.1 平流试验设计和数值计算方案

试验设计如图 8所示，有一气团初始时刻(t=0 s)密 度ρ分布满足

式中，，密度ρ₀=1，气团

半径宽R_x=25 km，高R_z=3 km，气团位于(x₀，z₀)=(-50 km，9 km)。无误差的情况下：平流运动5000 s后，气团位于山头之上，即(x₂，z₂)=(0 km，9 km)；平流运动10000 s后，气团位于山头的另一侧，与初始时刻的气团位置相对称，即(x₃，z₃)=(50 km，9 km)。图 8为平流试验的密度解析值分布，给定气团的平流速度为

式中，u₀=10 m/s，z₁=4 km，z₂=5 km。

质量通量方程的差分计算采用Arakawa-c网格(Arakawa et al，1977)和有限体积的一阶迎风方案。令式(16)中，对通量项采取分维算法

式(19)中通量计算公式为式中，(u，0)，u^-=Min(u，0)。 4.2 二维质量平流试验结果分析

由上述分析可知，在有地形情况下，不同坐标的坐标面起伏程度(亦称坡度)不同(图 1)，因而对质量输送的影响也不相同，可以通过比较输送过程中某时刻的密度场及其变化状况来说明不同坐标对质量输送的影响。根据Williamson等(1992)提出的计算全球积分误差的方法，本研究定义衡量误差公式为

式中，(ρ_i，k)num为数值计算结果，(ρ_i，k)ana为解析值，i、k为水平和垂直方向的格点下标。

下面分别比较在有地形和无地形情况下，t=0、5000、10000 s时4种坐标的密度分布，t=10000 s时的质量平流误差分布，以及积分各个时刻质量平流标准化误差演变。

图 9是一阶迎风方案下的4种坐标质量平流试验结果，其中图 9a₆、b₆是无地形数值模拟结果，可作为该数值计算方案的参考解。图 9a₁—a₆为不同坐标下对应于图 8的3个时刻的密度场模拟结果。与图 8的解析值比较可以看出，一阶迎风方案下无地形时质量气团随时间积分出现了水平拉伸，这主要是数值计算方案的耗散性决定的。COS坐标(Z_c=10 km)与无地形结果最接近，这是在临界衰减高度Z_c=10 km以下该坐标对地形平滑的结果；SLEVE2与COS坐标(Z_c=15 km)在t＝5000 s时气团形变相对于Gal.C.S与SLEVE1坐标都较小，只在气团的下部有一些地形影响特征；Gal.C.S与SLEVE1坐标在经过地形上方时气团形变较严重，并带有明显的地形影响特征，t＝10000 s时气团形状与无地形时相差甚远。图 9b₁—b₆为t＝10000 s的质量平流误差(即该时刻密度数值解与解析解之差)分布。4种坐标与无地形时在该时刻有不同的质量平流误差，且误差从Gal.C.S、SLEVE1、SLEVE2到COS坐标逐渐减小；无地形情况下绝对误差最大值为0.22，Gal.C.S、SLEVE1、SLEVE2坐标绝对误差最大值分别为0.65、0.61、0.38，而COS坐标(Z_c=10，15 km)与无地形最接近，误差最大值为0.32；特别是COS(Z_c=15 km)在较高层比SLEVE2坐标绝对误差小。

图 10是4种坐标和无地形情况下标准化误差的时间演变。图 10a为l₂误差演变，可以看出，4种坐标下平流输送的l₂标准化误差均随时间增大，在2000—7000 s误差迅速增大，这是因为2000—7000 s气团恰好通过地形上方，可以看出，地形对平流误差的影响较大；而7000 s以后误差变化不大。相比较而言，SLEVE2和COS坐标的误差是4种坐标中较小的，且COS坐标Z_c=10 km比Z_c=15 km时误差小，Z_c=10 km时误差最接近无地形情况。图 10b为误差演变l_∞情况，其误差分布与图 10a大致类似，它反映了最大误差(极值)随时间的演变。特别是COS坐标Z_c=10 km时误差与无地形时基本重合。不论是l₂还是l_∞，10000 s时Gal.C.S坐标的标准化误差比无地形大50%左右，而SLEVE2和COS坐标相对更接近无地形情况。

由试验结果和误差分析可以看出，平流误差与坐标面的起伏有关。SLEVE2、COS与Gal.C.S、SLEVE1坐标相比越往高层坐标面越平缓，因而可以更好地减小平流误差；并且，COS坐标临界衰减高度Z_c越小越接近无地形时的误差。5 结论和讨论

选取了4种典型的高度地形追随坐标，将其化为统一的坐标形式，分析其地形衰减系数b、坐标转换的雅可比项J_b以及坐标面上气压梯度力计算误差，然后分别进行了二维质量平流试验。通过理论分析和理想试验，可以得到以下几个初步结论：

(1)高度地形追随坐标可以转化为统一形式z=+bZ_S，将对坐标的研究简化为对地形衰减系数b的研究。

(2)SLEVE2坐标要满足坐标转换的J_b^-1>0条件，需要根据实践和理论合理地组合地形高度的各个尺度谱段Z_Si及与其对应的特征尺度h^*_i(i=1，2)。

(3)COS坐标的J_b^-1受到地形高度Z_S、衰减临界高度Z_c和参数n的约束。要满足坐标转换的J_b^-1>0条件：Z_S和Z_c确定时，n取值越大坐标面衰减效果越好，n取值有上界，且上界范围内n也不宜取值过大；Z_c和n确定时，地形高度也有上界，若取临界衰减高度Z_c为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坐标中临界衰减高度Z_c的选择很重要，Z_c越小越接近无地形的结果。

以上讨论的结果是基于理论分析和理想数值试验，各种高度地形追随坐标在真实大气数值预报模式如GRAPES模式(张人禾等，2008;陈德辉等，2008;杨学胜等，2008)中的表现还需要实际检验；同时，可以考虑结合各种坐标的优点设计新的高度地形追随坐标。