地形追随坐标自从由Phillips(1957)首次引入大气数值模式以来，由于其极大地简化了下边界条件的处理，因此得到了广泛的应用(Kasahara，1974; McNider et al，1981; Yamada，1983)。地形追随坐标分为气压(质量)地形追随坐标和高度地形追随坐标。两种坐标各有优缺点，从计算上来讲，气压地形追随坐标计算稳定性要好(Bénard，2003)，而高度地形追随坐标与时间无关，其尺度转换项随时间是常数，在半隐式时间积分中使用十分方便(Luo et al，2004)。在高分辨率可压缩非静力模式中为了较好地考虑垂直声波和重力波等快波，一般采用高度地形追随坐标。

但是地形追随坐标也带来了新的问题，气压梯度力的计算由一项变成了两个大项的小差，容易产生较大的计算误差，特别是在地形陡峭处这种误差变得尤其明显(Mahrer，1984；Dempsey et al，1998)，可产生强烈的虚假环流，从而影响模式的预报精度，甚至导致模式积分不稳定。另外，垂直坐标的转换所引入的计算误差不可避免，需要保证其尽量不影响大尺度动力特征(一般在高层)。为解决高度地形追随坐标下的气压梯度力计算误差问题，一般从两个方面入手，一是采用的垂直坐标使得地形对模式坐标面的影响随高度的增高尽量减小，二是采用高精度的气压梯度力计算方案。

中国气象局的全球和区域同化预报多尺度统一模式(GRAPES)所采用的Gal-Chen等(1975)提出的高度地形追随坐标，由于其形式简单，得到了广泛的应用，但该坐标衰减速率慢，地形的影响可一直延伸至模式高层(Leuenberger et al，2010; Klemp，2011)。Schör等(2002)提出了平缓地形追随坐标(Smooth Level Vertical coordinate，SLEVE)，该坐标将地形分为大小两种尺度，分别使用不同的衰减系数，使得小尺度地形的影响随高度增高很快地被移除，但该坐标对地形的衰减具有一定的尺度选择性。李超等(2012)将平缓地形追随坐标SLEVE坐标引入了GRAPES模式，对一次强降水天气过程进行了数值模拟，对比分析了和Gal-Chen坐标的影响，指出SLEVE坐标通过平滑坐标面上的地形作用，减小了虚假垂直速度进而减小虚假降水。Klemp(2011)提出了一种更为复杂的高度地形追随坐标，该坐标需要对每一层坐标面的高度进行平滑，从而滤去小尺度地形的影响，其效果在几个全球非静力模式的应用中得到了证实(Skamarock，2012)，下一代的WRF模式也将采用该坐标^①。

①Jimy Dudhia，个人通信。

为了减少高度地形追随坐标下气压梯度力计算误差问题，Mahrer(1984)提出将气压线性插值到与速度点相同的高度，然后计算水平气压梯度力，避免了地形修正项的计算，进而可以提高计算精度。Dempsey等(1998)指出在插值气压时使用二次多项式插值可以得到更高的计算精度。而在提出新的高度地形追随坐标的同时，Klemp(2011)也给出了同类较为简单的气压梯度力计算方案。

本研究将Klemp(2011)和Schör等(2002)提出的高度地形追随坐标分别引入GRAPES区域模式，通过理想实验，与现有Gal-Chen坐标进行比较，考察各种坐标下地形对高层大气的影响，并且评估不同坐标下气压梯度力的计算精度，为业务模式选择垂直坐标提供理论与实践基础。同时，为进一步减小气压梯度力的计算误差，引进了高精度的气压梯度力计算方案，以期提高GRAPES模式动力过程计算的精度与稳定性。

2 高度地形追随坐标介绍

考虑了3维科里奥利力与水平局地旋转的GRAPES模式运动方程在高度地形追随坐标下可写成

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

A()

(1)取A()=1- /z_T 时，A()随高度线性衰减，对应的坐标为Gal-Chen等(1975)提出的基本地形追随坐标，简称为Gal-Chen坐标，其中z_T为模式层顶高度。

(2)地形衰减廓线取为

点击浏览公式

(3)Klemp(2011)平滑坐标面的地形追随坐标可表示为 z=

点击浏览公式

h=h(x，y，)

点击浏览公式

点击浏览公式

点击浏览公式

A()

点击浏览公式

追随坐标，但它的最主要特点是不对地形进行尺度分离，而是对小地形进行水平滤波处理，随高度的增高逐层地消除小地形的影响，同时采用了混合地形衰减廓线A()。

图 1 为3种不同地形追随坐标所对应的坐标面()。SLEVE坐标的衰减特征高度取s= 8 km，Klemp坐标衰减特征高度z_H=14 km。可看到Gal-Chen坐标地形对模式面的影响可以一直延伸至模式顶层，SLEVE坐标虽有所改善，但仍不明显，Klemp坐标的改善最为显著，在高层模式坐标面趋于平直。

图 1 3种不同垂直坐标所对应的坐标面 (a. Gal-Chen坐标，b. SLEVE坐标，c. Klemp坐标;模式顶为zT=20 km，水平分辨率Δx= 500 m，垂直分辨率Δ= 500 m，垂直分40层) Fig. 1 Vertical coordinate surfaces for the (a) Gal-Chen coordinate，(b)SLEVE coordinate，and (c)Klemp coordinate implemented in the GRAPES model(The model top zT is 20 km， Δx=500 m，Δ= 500 m，with 40 vertical levels)

图选项

下面通过静止大气试验和重力波试验来评估3种坐标对气压梯度力计算精度以及重力波模拟能力的影响。以下试验均是基于GRAPES模式动力框架。 3.1 静止大气试验

静止大气试验的理想地形采用类似于Schör等(2002)的廓线

点击浏览公式

图 2为模式积分6 h后的位温垂直剖面。虽然低层有较强的稳定层结，Gal-Chen坐标中模式地形仍可在中高层产生明显影响，对位温场造成显著的扰动(图 2a)；SLEVE坐标地形同样在中高层有影响，但较Gal-Chen坐标有了明显的改善(图 2b)。效果最好的为Klemp坐标，对位温面产生的扰动最小，中高层位温面近于平直(图 2c)。

图 2 静止大气试验积分6 h后位温场的垂直剖面 (a. Gal-Chen坐标，b. SLEVE坐标，c.Klemp坐标;等值线间隔为1 K，地形未标出。地形为Y方向上不变化的带状地形) Fig. 2 Potential temperature as a function of height(contour interval 1 K)at t= 6 h from the resting-atmosphere simulation for the Gal-Chen coordinate(a)，SLEVE coordinate(b)，and Klemp coordinate(c)

图选项

静止大气试验通常可用于检验气压梯度力的计算精度。如果在理想状态下，模式没有误差，用水平均一的热力场初始化静止大气，模式在积分过程中大气应始终保持静止状态。但在实际数值模式中，由于计算误差的存在，模式大气会产生虚假运动，所以可以用垂直速度的大小来衡量计算误差的大小。图 3为静止大气试验中6 h积分时段内全场最大垂直速度随时间的变化。可以看出，Klemp坐标所产生的最大垂直速度始终小于1 m/s，而Gal-Chen坐标最大垂直速度甚至达到了4 m/s，SLEVE坐标介于两者之间。从以上理想试验可以看出，相比Gal-Chen和SLEVE坐标，Klemp坐标无论从模式地形对大气上层的影响，还是气压梯度力计算精度来讲，都具有最佳的效果。

图 3 静止大气试验3种不同坐标最大垂直速度的时间序列 Fig. 3 Time series of the maximum vertical velocity for the resting-atmosphere simulation for the Gal-Chen coordinate(thick solid line)，SLEVE coordinate(dashed line) and Klemp coordinate(thin solid line). The model top zT is 20 km

图选项

为了凸出垂直坐标对地形重力波造成的影响，在重力波试验中采用孤立圆形山作为理想地形，地形廓线取为

点击浏览公式

式中，山体最大高度h_m=500 m，山体宽度半径为a=3000 m。模式层顶高度为z_T=21 km，水平分辨率Δx=1000 m，垂直分辨率Δ= 350 m，垂直分60层。假定初始时刻大气满足静力平衡关系，过山气流风速(u)为8 m /s，地面气温取为290 K，浮力频率(N)为0.012 s^-1。可以发现无量纲量Na/u≈5＞1，根据Smith(1980)线性理论分析，这时气流经过山体时地形激发的重力波为静力(平衡)重力波。

重力波试验中，模式积分2 h后结果基本平衡，图 4给出了3种垂直坐标重力波试验积分3 h后的垂直速度高度剖面。与解析解相比，3种坐标都模拟出了静力重力波的垂直向上倾斜传播，但Gal-Chen坐标下的重力波表现出严重的扭曲和变形，垂直速度场有较大误差，Klemp坐标下对重力波的刻画与解析解最为接近，而SLEVE坐标介于两者之间。根据Smith(1980)对稳定层结下定常气流通过孤立山体时对气流扰动的线性理论分析，可以调整无量纲量Na/u≈1，从而产生非静力重力波，Klemp坐标下对非静力重力波的刻画同样优于Gal-Chen和SLEVE坐标(图略)。

图 4 积分3 h后3种不同垂直坐标的静力重力波试验垂直速度剖面 (a. Gal-Chen坐标，b. SLEVE坐标，c. Klemp坐标; 实线表示垂直向上的速度，虚线表示垂直向下的速度，等值线间隔为0.06 m/s) Fig. 4 Vertical cross section of vertical velocity for the hydrostatic gravity wave test after 3 hours (a)Gal-Chen coordinate，(b)SLEVE coordinate，and (c)Klemp coordinate. Contour intervals are 0.06 m/s. The downward(upward)vertical velocity is in the dashed(solid)line

图选项

本节将针对气压梯度力计算方案做进一步改进，以提高其计算精度。

在模式计算中，x方向的气压梯度力从z坐标至高度地形追随坐标的转换关系为

点击浏览公式

可见气压梯度力由一项变为了两项，右边第一项为模式坐标上的气压梯度力，第二项为地形追随坐标垂直变化的补偿项。

使用Gal-Chen坐标，经过众多不同水平和垂直分辨率的试验，发现气压梯度力计算误差所导致虚假环流的强度对水平分辨率不是十分敏感，而强烈地依赖于模式的垂直分辨率，即随着垂直格距Δ的减小而减小(图 5)，说明气压梯度力的计算误差主要来源于式(21)右边补偿项的影响。为了避免该项的计算，参考Mahrer(1984)方法，Klemp(2011)提出了一种简化的Mahrer方法，该方法适用于地形相对平缓或模式垂直分层相对稀疏(图 6a)情况下气压梯度力的计算，不需要判定与水平速度点具有相同高度的点位于哪两个坐标面之间，可直接利用邻近的两层气压进行插值，易于编程，且不需要增加额外的计算量。但是，在陡峭地形下，对与速度点相同高度的水平气压进行插值时，往往需要超过上下相邻两层的气压，该简化的方法可能会引起更大的计算误差。因此，与Klemp(2011)中简化的Mahrer方法不同，本研究采用更为一般的Mahrer方法(图 6b)，需判断与水平速度点具有相同高度的 气压点位于哪两个坐标面之间，然后根据相邻的两层气压进行插值，这种方法更适合相对陡峭的地形和一般的垂直分层，更具有普遍性。Dempsey等(1998)指出一般的Mahrer方法计算气压梯度力的精度最高。

图 5 不同垂直分层的静止大气试验最大垂直速度的时间序列(模式顶为zT=20 km，水平分辨率Δx=500 m) Fig. 5 Time series of the maximum vertical velocity for the resting-atmosphere simulation for the different number of vertical level: 40 levels(thick solid line)，50 levels(dashed line)，60 levels(thin solid line).The model top zT is 20 km，Δx=500 m

图选项

图 6 倾斜坐标面下计算水平气压梯度力(a)简化的Mahrer方法和(b)一般情况的Mahrer方法示意图(图b中插值点由箭头指示) Fig. 6 Schematical charts for computing the horizontal pressure gradient force at sloping coordinate surfaces using a simplified version of Mahrer's approach(a)，and a generalized Mahrer's approach(b)(for further details see the text)

图选项

在水平Arakawa-C网格下，高度地形追随坐标中一般的Mahrer气压梯度力计算方法为

点击浏览公式

式中，Π^*是在速度点相同高度面上的气压值，由相邻的两层气压插值得到。当Π^*的高度位于地面以下或模式顶层以上时，则由最底3层模式层或最高3层模式层上的气压值做2次多项式的拉格朗日插值得到。与式(21)不同，式(22)避免了垂直变化补偿项的计算。

图 3表明Klemp坐标本身已对气压梯度力的计算精度有了很大的改进，为了更能说明问题，仅用Gal-Chen坐标检验该方案对气压梯度力的改进作用，与上节类似，设计静止大气试验来测试新气压梯度力计算方案的精度。Mahrer(1984)指出，当模式坐标面的坡度大于模式格距比Δz/Δx时，水平气压梯度力的计算误差显著增大。通常模式的水平分辨率都是远大于垂直分辨率的，因此为了凸出一般的Mahrer气压梯度力计算方案的精度，水平分辨率取Δx= 500 m，采用非均匀分层，逐渐增加垂直层数，分别对垂直51层、71层和81层进行了试验。图 7为Gal-Chen垂直坐标下的新老气压梯度力计算精度的比较。从中可以看出，新的气压梯度力计算方案对气压梯度力的精度有较大的改进。

图 7 Gal-Chen坐标下不同垂直分层的静止大气试验最大垂直速度的时间序列 (实线为采用式(21)直接计算气压梯度力，虚线为采用式(22)一般的Mahrer方法计算气压梯度力; a. 垂直51层，b. 垂直71层，c. 垂直81层；垂直方向采用非均匀分层；模式顶为zT=20 km，水平分辨率Δx=500 m，水平计算区域为60 km×25 km) Fig. 7 Time series of the maximum vertical velocity for the resting-atmosphere simulation with vertically stretched grids using the Gal-Chen coordinate for the different numbers of vertical level:(a)51 levels，(b)71 levels，and (c)81 levels. (Solid lines are ones that use the original formulation for the horizontal pressure gradient force，dashed lines the generalized Mahrer's approach. The model top zT is 20 km，and Δx=500 m. The computational domain is 60 km×25 km)

图选项

为简单起见，使用Gal-Chen坐标和Klemp坐标分别对实际天气个例进行模拟，通过分析物理量场分布比较两种坐标的优劣。模式水平分辨率为0.1°，垂直51层，起报时间为2014年5月20日08时(北京时)，进行24 h预报。图 8a和b分别为Gal-Chen坐标和Klemp坐标下沿35°N的位温垂直剖面，由图可见，Klemp坐标的位温场等值线在高层明显平滑，等值线上由地形引起的小尺度扰动有很大程度的减小(图 8b)。而Gal-Chen坐标的位温场在高层依然受地形的影响(图 8a)，与前面理想试验的结果基本一致。图 8c和d分别给出了Gal-Chen坐标和Klemp坐标试验模式第20层的24 h位温预报场。从图可以明显看出，Gal-Chen坐标产生很强的位温扰动(图 8c)，且扰动与小尺度地形密切相关，而使用Klemp坐标的位温场就平滑的多(图 8d)，可见其对地形影响的衰减更加明显。

图 8 24 h预报的实际个例的沿35°N位温垂直剖面 (a、b.等值线间隔为20 K; a. Gal-Chen坐标，b. Klemp坐标; c、d. 模式第20层的位温水平分布(等值线间隔为5 K); c. Gal-Chen坐标，d.Klemp坐标) Fig. 8 24 h forecast results of a real case simulation for (a)Gal-Chen coordinate，and (b)Klemp coordinate. The above panel shows that potential temperature along the cross section at 35°N(Contour intervals are 20 K). The bottom panel shows the potential temerature horizontal fields at the 20th model level for(c)Gal-Chen coordinate，and (d)Klemp coordinate(Contour intervals are 5 K)

图选项

将一种新的高度地形追随坐标(Klemp坐标)和高精度的气压梯度力方案引入了GRAPES区域模式，通过静止大气试验、重力波试验及实际个例模拟对结果进行了评估，并与Gal-Chen 坐标和SLEVE坐标进行了比较，得到如下结论：

(1) Klemp坐标相比Gal-Chen和SLEVE坐标，对气压梯度力计算精度有较大的改进，静止大气试验中Klemp坐标的最大垂直速度始终小于其他两种坐标，积分过程中模式高层的位温面基本平直，没有出现较大的虚假扰动。

(2) 相比其他两种坐标，Klemp坐标对静力重力波的刻画最佳，结果与解析解接近。Gal-Chen和SLEVE坐标对重力波的模拟出现了不同程度的扭曲和变形。

(3) 与Klemp(2011)简化的Mahrer方案不同，将更为一般的Mahrer方案引入GRAPES用来计算水平气压梯度力，试验证明该方案对气压梯度力的计算精度有相当的提高。该方法更适合相对陡峭的地形和更一般的垂直分层，但需要判断与速度点具有相同高度的气压点位于哪两个坐标面之间。

与SLEVE坐标一样，Klemp坐标在具体的应用过程中，需要对各种参数，如β_k和z_H进行估计，以确定最优的数值。由于Klemp坐标在每一模式层上都对地形进行了一定程度的平滑，实际上在每一模式层上具有不同的h(x，y，)，因此在实际应用于GRAPES模式当中，需要将地形坡度变量由2维变量变为3维变量，由此需要对GRAPES程序进行较大的修改。

简化的Mahrer计算方案适合中等坡度的地形或模式垂直分层不够密集的情况，如果地形非常陡峭或者模式低层垂直分层密集，简化的Mahrer方案的精度有所下降，因此采用了更一般的Mahrer方法，需要对与速度点具有相同高度的气压点位于哪两个坐标面之间进行判断。同时，在插值与速度点在相同高度上的气压时，采用的是线性插值，为了提高计算精度可以采用2阶插值方案，因此在地形非常陡峭处气压梯度力计算精度的问题值得进一步研究和探讨。

以上的结果大部分基于理想试验，实际个例的分析局限性较大，在真实大气数值模拟中尤其在应用于GRAPES区域模式后的预报效果需要做进一步的检验。