中国气象学会主办。
文章信息
- 张旭, 黄伟, 陈葆德. 2015.
- ZHANG Xu, HUANG Wei, CHEN Baode. 2015.
- 一种新型高度地形追随坐标在GRAPES区域模式中的实现:理想试验与比较研究
- Implementation of the Klemp height-based terrain-following coordinate in the GRAPES regional model: Idealized tests and inter-comparison
- 气象学报, 73(2): 331-340
- Acta Meteorologica Sinica, 73(2): 331-340.
- http://dx.doi.org/10.11676/qxxb2015.014
-
文章历史
- 收稿日期:2014-04-14
- 改回日期:2014-10-10
2. 中国气象局上海台风研究所, 上海, 200030
2. Shanghai Typhoon Institute of China Meteorological Administration, Shanghai 200030, China
地形追随坐标自从由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模式运动方程在高度地形追随坐标下可写成
式中,φ0为旋转经纬度坐标系统下模式北极点的地理坐标,而 是垂直坐标的转换项,且满足 与单调可逆条件高度地形追随坐标可统一写成式中,A()为地形特征的衰减率,zs为地形高度。取不同的地形衰减廓线则对应不同的高度地形追随坐标。(1)取A()=1- /zT 时,A()随高度线性衰减,对应的坐标为Gal-Chen等(1975)提出的基本地形追随坐标,简称为Gal-Chen坐标,其中zT为模式层顶高度。
(2)地形衰减廓线取为
对应的坐标为Schör等(2002)提出的平缓地形追随坐标。该坐标需要将地形分为大尺度地形和小尺度地形,对不同的地形特征采用不同的衰减特征高度(s),在此高度内地形特征大约衰减1/e,s的取值应根据实际地形特征进行确定。以下试验中不对地形进行尺度分离,即只取单一的衰减特征高度。(3)Klemp(2011)平滑坐标面的地形追随坐标可表示为 z=
式中,h=h(x,y,)为实际高度,且h(x,y,0)=zs(x,y)为实际地形高度,该坐标要求在不同的坐标高度()上对实际高度进行不同程度的平滑。平滑方案可取为 实际应用中,用下式计算 式中,n为平滑迭代次数,一般的取值范围为20—30;且 式中,hm为最大地形高度。对于混合地形衰减廓线A()有 式中,zH为衰减特征高度。Klemp坐标与SLEVE坐标相同,本质上也是一种混合高度地形追随坐标,但它的最主要特点是不对地形进行尺度分离,而是对小地形进行水平滤波处理,随高度的增高逐层地消除小地形的影响,同时采用了混合地形衰减廓线A()。
图 1 为3种不同地形追随坐标所对应的坐标面()。SLEVE坐标的衰减特征高度取s= 8 km,Klemp坐标衰减特征高度zH=14 km。可看到Gal-Chen坐标地形对模式面的影响可以一直延伸至模式顶层,SLEVE坐标虽有所改善,但仍不明显,Klemp坐标的改善最为显著,在高层模式坐标面趋于平直。
3 理想试验下面通过静止大气试验和重力波试验来评估3种坐标对气压梯度力计算精度以及重力波模拟能力的影响。以下试验均是基于GRAPES模式动力框架。 3.1 静止大气试验
静止大气试验的理想地形采用类似于Schör等(2002)的廓线
式中,a = 5 km,λ= 4 km,地形最大高度hm=1000 m。为了在低层包含较为复杂静止大气的垂直气压梯度结构,在一致静力稳定的大气层结(浮力频率取N=0.01 s-1)结构中,在2—3 km加入一逆温层(N=0.02 s-1);初始气压场由静力平衡关系给定。参考大气廓线为绝热大气θ=288 K。大气的初始状态为u=0,v=0,w=w=0。模式顶为zT=20 km,水平分辨率Δx=500 m,垂直分辨率Δ= 500 m,垂直分40层。注意理想地形的廓线式(19)只给出了二维形式,在试验中采用的是三维的带状地形,地形在y方向上没有变化。
图 2为模式积分6 h后的位温垂直剖面。虽然低层有较强的稳定层结,Gal-Chen坐标中模式地形仍可在中高层产生明显影响,对位温场造成显著的扰动(图 2a);SLEVE坐标地形同样在中高层有影响,但较Gal-Chen坐标有了明显的改善(图 2b)。效果最好的为Klemp坐标,对位温面产生的扰动最小,中高层位温面近于平直(图 2c)。
静止大气试验通常可用于检验气压梯度力的计算精度。如果在理想状态下,模式没有误差,用水平均一的热力场初始化静止大气,模式在积分过程中大气应始终保持静止状态。但在实际数值模式中,由于计算误差的存在,模式大气会产生虚假运动,所以可以用垂直速度的大小来衡量计算误差的大小。图 3为静止大气试验中6 h积分时段内全场最大垂直速度随时间的变化。可以看出,Klemp坐标所产生的最大垂直速度始终小于1 m/s,而Gal-Chen坐标最大垂直速度甚至达到了4 m/s,SLEVE坐标介于两者之间。从以上理想试验可以看出,相比Gal-Chen和SLEVE坐标,Klemp坐标无论从模式地形对大气上层的影响,还是气压梯度力计算精度来讲,都具有最佳的效果。
3.2 重力波试验为了凸出垂直坐标对地形重力波造成的影响,在重力波试验中采用孤立圆形山作为理想地形,地形廓线取为
式中,山体最大高度hm=500 m,山体宽度半径为a=3000 m。模式层顶高度为zT=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 气压梯度力计算方案的改进本节将针对气压梯度力计算方案做进一步改进,以提高其计算精度。
在模式计算中,x方向的气压梯度力从z坐标至高度地形追随坐标的转换关系为
可见气压梯度力由一项变为了两项,右边第一项为模式坐标上的气压梯度力,第二项为地形追随坐标垂直变化的补偿项。
使用Gal-Chen坐标,经过众多不同水平和垂直分辨率的试验,发现气压梯度力计算误差所导致虚假环流的强度对水平分辨率不是十分敏感,而强烈地依赖于模式的垂直分辨率,即随着垂直格距Δ的减小而减小(图 5),说明气压梯度力的计算误差主要来源于式(21)右边补偿项的影响。为了避免该项的计算,参考Mahrer(1984)方法,Klemp(2011)提出了一种简化的Mahrer方法,该方法适用于地形相对平缓或模式垂直分层相对稀疏(图 6a)情况下气压梯度力的计算,不需要判定与水平速度点具有相同高度的点位于哪两个坐标面之间,可直接利用邻近的两层气压进行插值,易于编程,且不需要增加额外的计算量。但是,在陡峭地形下,对与速度点相同高度的水平气压进行插值时,往往需要超过上下相邻两层的气压,该简化的方法可能会引起更大的计算误差。因此,与Klemp(2011)中简化的Mahrer方法不同,本研究采用更为一般的Mahrer方法(图 6b),需判断与水平速度点具有相同高度的 气压点位于哪两个坐标面之间,然后根据相邻的两层气压进行插值,这种方法更适合相对陡峭的地形和一般的垂直分层,更具有普遍性。Dempsey等(1998)指出一般的Mahrer方法计算气压梯度力的精度最高。
在水平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垂直坐标下的新老气压梯度力计算精度的比较。从中可以看出,新的气压梯度力计算方案对气压梯度力的精度有较大的改进。
5 不同垂直坐标的实际个例的模拟分析为简单起见,使用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),可见其对地形影响的衰减更加明显。
6 总结与讨论将一种新的高度地形追随坐标(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和zH进行估计,以确定最优的数值。由于Klemp坐标在每一模式层上都对地形进行了一定程度的平滑,实际上在每一模式层上具有不同的h(x,y,),因此在实际应用于GRAPES模式当中,需要将地形坡度变量由2维变量变为3维变量,由此需要对GRAPES程序进行较大的修改。
简化的Mahrer计算方案适合中等坡度的地形或模式垂直分层不够密集的情况,如果地形非常陡峭或者模式低层垂直分层密集,简化的Mahrer方案的精度有所下降,因此采用了更一般的Mahrer方法,需要对与速度点具有相同高度的气压点位于哪两个坐标面之间进行判断。同时,在插值与速度点在相同高度上的气压时,采用的是线性插值,为了提高计算精度可以采用2阶插值方案,因此在地形非常陡峭处气压梯度力计算精度的问题值得进一步研究和探讨。
以上的结果大部分基于理想试验,实际个例的分析局限性较大,在真实大气数值模拟中尤其在应用于GRAPES区域模式后的预报效果需要做进一步的检验。
李超, 陈德辉, 李兴良. 2012. 大气数值预报模式中高度地形追随坐标的设计研究: 理论分析与理想试验. 气象学报, 70(6): 1247-1259. Li C, Chen D H, Li X L. 2012. A design of height-based terrain-following coordinates in the atmospheric numerical: Theoretical analysis and idealized tests. Acta Meteor Sinica, 70(6):1247-1259 (in Chinese) |
BénardP. 2003. Stability of semi-implicit and iterative centered-implicit time discretizations for various equation systems used in NWP. Mon Wea Rev, 131(10): 2479-2491 |
Dempsey D, Davis C. 1998. Error analyses and tests of pressure-gradient force schemes in a nonhydrostatic, mesoscale model//12th Conferenceon Numerical WeatherPrediction.Phoenix, AZ:American Meteorological Society, 236-239 |
Gal-Chen T, Jsomerville 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 |
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 |
Leuenberger D, Koller M, Futher O, et al. 2010. A generalization of the SLEVE vertical coordinate. Mon Wea Rev, 138(9): 3683-3689 |
Luo H X, Bewley T R. 2004. On the contravariant form of the Navier-Stokes equations in time-dependent curvilinear coordinate systems. J Comput Phys, 199(1): 355-375 |
Mahrer Y. 1984. An improved numerical approximation of the horizontal gradients in a terrain-following coordinate system. Mon Wea Rev, 112(5): 918-922 |
McNider R T, Pielke R A. 1981. Diurnal boundary-layer development over sloping terrain. J Atmos Sci, 38(10): 2198-2212 |
Phillips N A. 1957. A coordinate system having some special advantages for numerical forecasting. J Meteor, 14(2): 184-185 |
Schör C, LeuenbergerD, FuhrerO, et al. 2002. A new terrain-following vertical coordinate formulation for atmospheric prediction models. Mon Wea Rev, 130(10): 2459-2480 |
Skamarock W C, Klemp J B, DudaM G, et al. 2012. A multi scale nonhydrostatic atmospheric model using centroidal Voronoi tesselations and C-grid staggering. Mon Wea Rev, 140(9): 3090-3105 |
Smith R B. 1980. Linear theory of stratified hydrostatic flow past an isolated mountain. Tellus, 32(4): 348-364 |
Yamada T. 1983. Simulations of nocturnal drainages flows by a q2l turbulence closure model. J Atmos Sci, 40(1): 91-106 |