1 引言

当层结稳定的气流翻越山脉时可能激发出地形重力波 (Durran, et al, 1982)。地形重力波向上传播的过程中，由于密度减小而加深其振幅, 从而导致波的破碎 (Eliassen, et al, 1961)，此外，当相对于波动尺度地形高度足够高时，重力波进入到非线性过程中，气流绕行在山后产生扰动也会导致波破碎，重力波破碎会移除波的能量引起大尺度流场的耗散作用，这种耗散作用称为地形重力波拖曳。地形重力波的水平尺度约为数千米至百千米，而时间尺度在数小时至一天，对于现在的多数数值模式而言，尚不能完整地解析这一现象。

Boer等 (1984)第一次在大尺度模式中考虑了次网格地形重力波拖曳 (GWD) 对低层大气环流的影响。Palmer等 (1986)基于Lindzen (1981)的饱和假设理论，发展了GWD参数化方案，该方案主要考虑平流层低层波破碎对重力波拖曳的影响，发现参数化后对流层西风急流偏强的问题得以减弱。McFarlane (1987)在气候模式中加入了GWD参数化方案，发现地形激发的重力波破碎造成的动量下沉对于对流层和平流层低层之间的气流的结构起着决定性作用。Kim等 (1995)在二维滞弹性非静力的中尺度模式中加入GWD参数化方案 (简称KA95方案)，结果表明采用GWD参数化，大气低层重力波的破碎可使得波能在低层被捕捉，并通过非静力波的共振使得拖曳力增强。随后在气候模式 (Kim, 1996) 和区域模式 (Alpert, et al, 1996) 中运用和发展了这种参数化方案。Lott等 (1997)提出的GWD参数化方案中，认为当地形有效高度足够高时，该方案可以有效地处理低层的阻塞气流，其预报结果更加符合实际观测。Smith等 (2006)用7个气流过山例子探讨了模式水平分辨率对地形引起的重力波拖曳的敏感度，结果表明在多数个例中，即使水平分辨率提高至4 km以上，仍无法完整描述模式中GWD的影响。

高守亭等 (2003)对McFarlane (1987)的方案进行了改进，得到了较完整的描述定常重力波破碎对纬向平均流拖曳作用的参数化方案，该方案不仅考虑了重力波破碎后尚存的动量通量对纬向平均流的拖曳作用，还考虑了在重力波破碎区耗散作用造成的动量损失对纬向平均气流的影响。徐国强等 (2010)在GRAPES全球模式中引入欧洲中期天气预报中心 (ECMWF) 的GWD参数化方案并进行了敏感性试验，结果表明，引进GWD过程后使预报的流场更接近于大气真实状态，从而提高了降水预报准确率。钟水新等 (2014)将KA95方案引入到华南中尺度模式GRAPES中，改善了模式对冬、春季节低层南风偏强和地面气温过高的现象。

基于以上论述可以发现，关于GWD参数化无论是对区域数值天气预报还是气候模拟均具有非常重要的作用。青藏高原是世界上海拔最高、地形最复杂的高原，青藏高原约占中国土地总面积的四分之一，平均海拔超过4000 m。研究表明，在夏季青藏高原的作用主要表现在对气流的绕流作用上，在冬季青藏高原的作用在爬流和绕流两个方面均很重要 (丁一汇，2005)。爬流会产生地形重力波拖曳，而绕流则会形成阻塞拖曳。但目前对青藏高原大地形激发的重力波拖曳的相关问题仍不十分清楚，存在以下需要进一步研究的问题：(1) 青藏高原大地形激发的重力波拖曳力的分布；(2) 爬坡和绕流高度大致范围是怎样的；(3) 青藏高原地形重力波拖曳对天气的影响。本研究首先在GRAPES区域中尺度数值预报系统 (GRAPES_Meso) 模式引入次网格地形重力波拖曳参数化方案，对以上问题进行了初步分析，希望对青藏高原地形重力波拖曳有进一步的认识。

2 模式和方案介绍 2.1 GRAPES_Meso模式

GRAPES为全球与有限区域通用、静力与非静力可选的多尺度数值预报模式。该系统垂直方向采用地形追随高度坐标，水平方向为球面坐标。在空间离散化时，模式水平方向取Arakwawa-C跳点经纬网格，垂直方向选择Charney-philips变量配置。时间积分使用半隐式半拉格朗日方案 (陈德辉等，2008；薛纪善等，2008)。模式开发建立了插拔式可移植的区域中尺度物理过程软件包，物理过程的完整性与国际上先进模式相当 (徐国强等，2008)。GRAPES_Meso为GRAPES区域中尺度数值预报系统，其核心部分是GRAPES预报模式系统包括动力框架以及经过优化和改进的物理过程参数化方案、三维变分资料同化系统、模式标准初始化系统。在文中个例试验和批量试验中，物理过程方案选择如表 1。

2.2 次网格地形重力波拖曳参数化方案

目前，GRAPES_Meso模式尚未考虑地形重力波拖曳的影响，因此首先在模式中引入了WRF中的次网格地形重力波拖曳参数化方案，该方案由Hong等 (2008)所发展，除了处理了传统重力波在高层的破碎 (Alpert, et al, 1988) 问题，也参考了KA95参数化方法，强化处理了因地形非对称性和尖锐性导致的低层波破碎的问题。此外，还结合了Lott等 (1997)提出的低层气流阻塞拖曳参数化方法 (简称LM97方案)。方案中主要流程描述如下：

山地的无量纲高度 (H_n)

(1)

式中，H为山脉的最大高度，U为水平风速，N为Brunt-vaisala频率。

当H_n较小时，气流爬过山脉并且激发地形重力波，这时主要考虑地形重力波拖曳。根据KA95方案，参考层 (即重力波产生的高度) 上产生的重力波拖曳

(2)

(3)

(4)

式中，ρ为密度，U为水平风速，Δx为模式水平分辨率，下标0代表参考层，E为参考层上拖曳力的增强因子，用来增加受到低层波破碎或是抑制作用导致的拖曳力，它会受到地形非对称性 (O_A) 以及弗劳德数 (Fr) 的影响，Fr_c为临界弗劳德数。m′代表网格内山的数量，受到地形非对称 (O_A) 以及地形有效长度 (L_x) 的影响。G′为一个用来断定气流阻塞与否的渐进函数 (Pierrehumbert，1996)，G_max为无量纲饱和通量 (一般设为1)。O_C代表地形尖锐性，C_E和C_g为常数。

得到参考层上的重力波拖曳之后，就可以计算拖曳力的垂直分布。重力波如果没有满足饱和假设，能量则可以在不损失的情况下向上传播，因此，使上层的重力波拖曳等于下层的重力波拖曳 (τ_i=τ_i+1)，从而来计算上层波的振幅 (h_di)，并利用计算的该振幅求取局地里查逊数 (Ri_m)

(6)

(7)

(8)

根据Lindzen (1981)饱和假设，当Ri_m < Ri_c(=0.25) 时，波会发生破碎。因此当没有达到饱和条件时，重力波在该层不会破碎，该层的拖曳力 (τ) 维持不变；反之，重力波破碎，部分能量损失，直到回到饱和临界点 (Ri_m=0.25)。因此可将Ri_m=0.25代入式 (9) 中计算临界振幅h_c

(9)

利用式 (9) 求得的振幅来计算该层的剩余重力波拖曳力

(10)

然后重复计算各层拖曳力分布直到τ_i=0或者到达模式层顶。

当H_n较大时，气流垂直运动受到限制，一部分气流可以爬坡越过山脉，一部分则阻塞形成绕流，产生阻塞拖曳力。阻塞层高度计算公式为

(11)

与爬坡产生的重力波拖曳不同，阻塞拖曳主要在阻塞层高度，因此，该阻塞拖曳力积分是从Z=0至Z=Z_b。

(12)

(13)

式中，l(z) 为与风向和阻塞高度有关的地形有效高度 (Alpert, 2004)，μ为山的高度取地形标准差，φ为低层风风向与地形主轴的夹角，参数a≈μ/σ，b≈a/γ，σ代表地形各向异性OD，γ代表地形坡度，C_d是拖曳系数，为常数。因此，单位面积的阻塞拖曳力为

(14)

(15)

在该方案中，将重力波拖曳和阻塞拖曳加起来，代入模式物理方程中。

2.3 模拟试验设计

通过是、否引入次网格地形重力波拖曳参数化方案进行两组对比试验，一组考虑参数化方案为敏感性试验 (GWD)，另一组不考虑参数化方案为控制试验 (CTRL)。两组试验采用的其他物理过程均相同，采取方案如表 1，模式模拟范围为 (22°—42°N，79°—119°E)，水平格距为0.1°，网格点数为401×201，垂直方向为不等距49层，时间分辨率60 s。通过GRAPES_Mesov4.0.0.0模式进行了个例对比分析和批量数值试验。文中选取的试验个例为2013年7月8—10日四川盆地一次降水过程。批量试验选取时段为2013年8月1—31日。

3 青藏高原地形重力波拖曳相关参数初步分析 3.1 个例天气形势分析

以2013年7月8日00时—10日00时 (世界时，下同) 的一次四川盆地强降水过程作为试验个例。这次降水过程覆盖了四川盆地的西北部地区 (图 6a、b)。这次强降水过程具有持续时间长、降水强度大、强降水落区主要位于盆地向高原过渡区的特点。从降水发生的前24 h对流层中层500 hPa的环流形势 (图略) 可以发现，在亚洲与欧洲大陆高纬度区域呈现东高西低的环流形势，西风气流位于中纬度地区，副热带高压北移至30°N，脊点位于120°E附近。在降水区西侧，与副热带高压同一纬度有伊朗高压相对应，这样在两个高压之间形成了十分宽广的槽，不断有短波槽从青藏高原移入到四川盆地。由于副热带高压的稳定维持，使得天气系统稳定少动或移动缓慢，槽前的偏南气流将孟加拉湾、南海水汽源源不断输送至盆地地区，为这次暴雨过程提供了充沛的水汽条件。

3.2 地形重力波拖曳力分布

当气流被山脉阻挡时，一部分气流可以爬过山脉，一部分则形成绕流。那么，对于这一次天气过程，青藏高原形成的重力波拖曳和阻塞拖曳是如何分布的呢？

图 1a、c分别为GWD试验模拟24 h后沿30°和32°N重力波拖曳和绕流阻塞拖曳垂直分布，其中暖色系阴影代表重力波拖曳，冷色系阴影代表阻塞拖曳，纵坐标为模式层层数；图 1b、d分别是对应的30°和32°N的地形剖面 (灰色阴影) 和模式层高度 (等值线)，其中红色等值线为模式层第5层，绿色等值线为第10层。可以看出，阻塞拖曳主要位于模式的第1—5层 (第5层距地表高度约为417 m)，大小为 (-1—-5)×10^-4N/m²；重力波拖曳主要位于第5—10层 (第10层距地表高度约为1440 m)，大小为 (-1—-6)×10^-3N/m²。由图 1a、b还可看出，在30°N，阻塞拖曳高度较高的地区，其对应的重力波拖曳高度也较高；阻塞拖曳较小的地区，其对应的气流爬坡激发的重力波拖曳较大。

图 2为GWD试验模拟24 h后模式层第3层 (距地表高度约182 m) 和第5层 (距地表高度约417 m) 重力波拖曳及阻塞拖曳的水平分布。从图 2a、b可以看出，第3层以阻塞拖曳为主，阻塞拖曳主要位于青藏高原边缘地区，阻塞拖曳大值区沿喜马拉雅山脉走向和青藏高原东坡，与地形标准差大的区域相对应，即地形起伏程度大的山脉对气流的阻挡作用更显著，从而绕流造成的阻塞拖曳也更显著。阻塞拖曳最大值位于高原南部30°N附近，达-7×10^-3 N/m²。模式层第5层以重力波拖曳为主，重力波拖曳主要位于青藏高原东部地区和云贵高原的北部边缘，最大值达-7×10^-3 N/m²，重力波拖曳较阻塞拖曳大一个量级。

3.3 气流爬坡和绕流

弗劳德数 (Fr) 是描述过山气流特征的一个重要参数, 其定义为Fr=U/Nh，其中U为气流平均速度，N为气层静力稳定度，h为山高。气流过山的浅水理论中，Fr决定过山气流的特征：Fr较小时，气流以绕流为主；Fr增大时，气流开始爬坡，即当大气层结稳定度一定时，风速越小、山越高，气流越容易绕流；风速越大、山越矮，气流越容易爬坡。一般情况下，当Fr>1时，气流就有了过山能力，在文中引入的参数化方案中，临界弗劳德数 (Fr_c) 取为1。Fr表明了气流过山的能力，对于高大的地形，气流爬升到一定高度之后就不能继续向上爬升从而产生绕流，这个高度就是绕流高度，也称为阻塞高度。

为对青藏高原的气流爬坡和绕流有更直观的认识，计算了弗劳德数和阻塞高度 (图 3)，图 3a为GWD试验模拟24 h后的Fr分布，在地势平坦的地区，如塔里木盆地、内蒙古高原、四川盆地和东部的平原地区，Fr均小于1，而Fr大于1的区域与地形标准差大的地区相对应，最大值在喜马拉雅山脉和高原东部边缘地区。图 3b为GWD试验模拟24 h后的绕流高度，绕流高度用模式层层数表示。对比图 3a、b可以看出，Fr大的地区即气流爬山能力强的地区，绕流高度也相应的比较高，沿喜马拉雅山脉和青藏高原东部边缘绕流高度最高，可达第7层 (距地表高度约为757 m)，青藏高原西部绕流高度在第3层左右 (距地表高度约为182 m)，中部则在第2层 (距地表高度约为107 m)。

3.4 重力波破碎

根据重力波破碎理论，在重力波向上传播的过程中，会因为密度的减小而加深其振幅，引起不稳定导致波的破碎，即高层重力波破碎；而低层重力波破碎是由于山峰表面和下游形成的山脉波的共振导致的。此外，当地形高度很高时，气流无法顺利爬升过地形时，会产生绕流，山后产生涡流，也会导致波的破碎。在重力波破碎的区域存在驻点即风速 (u) 快速减小接近0，等熵线也有剧烈变化，位温显著地减小。

图 4a、b、c分别为模拟24 h后沿32°N CTRL和GWD试验的位温和纬向风速 (u) 及两者纬向风速差的垂直分布。对比图 4a、b发现，采用次网格地形重力波拖曳参数化方案后，在中低层，83°—86°E山脉上空500—350 hPa位温变化比较剧烈，且风速 (u) 有明显的减小，表明可能存在重力波破碎现象。在104°—106°E, 800—500 hPa位温急剧减小，等位温线很陡，同时风速 (u) 也显著减小，尤其在700—550 hPa接近于0，表明在此下坡区域可能存在低层重力波破碎。在111°E附近，850—600 hPa等位温线也变陡，位温和风速 (u) 减小，说明在这一背风波位置同样可能发生重力波破碎。在高层 (超过300 hPa)，次网格地形重力波拖曳参数化之后，位温和风速 (u) 不像中低层那样变化剧烈，在93°和100°E的高层150—70 hPa位温和风速 (u) 减小，且在70 hPa附近u接近于0，表明此处可能存在高层重力波破碎。

由图 4c可以看出，两个试验的差值在垂直方向上呈正—负—正—负的波列形式，表明采用次网格地形重力波拖曳参数化方案之后，初步模拟出了重力波在垂直方向向上传播的特征，而且波列上传至70 hPa附近时消失，特别是在93°和100°E与图 4b中风速 (u) 和位温剧烈减小的区域相对应，表明在这一高度可能存在重力波破碎，与上述分析结论相吻合。

关于次网格地形重力波破碎问题这里只是进行了探索性地分析和研究，后续尚需更深入的研究和更多的观测事实来验证。

3.5 次网格地形重力波拖曳对风场和降水的影响 3.5.1 对风场的影响

次网格地形重力波拖曳和阻塞拖曳在模式物理过程中有引起大气水平运动场的倾向。为了定量评估采用次网格地形重力波拖曳参数化方案对于该次个例的风场模拟的改进效果，进行0—24和24—48 h风场检验对比。检验的区域为模拟区域。文中主要针对青藏高原及下游地区进行研究，因此选取西藏和四川的9个高空站进行对比分析 (表 2)。

图 5为次网格地形重力波拖曳试验和CTRL试验模拟24 h后风场与各高空站探空风的均方根误差，从图中可以看出采用次网格地形重力波拖曳参数化方案之后，拉萨站 (图 5a) 在中低层均方根误差都显著减小；那曲站 (图 5b) 各高度上均方根误差都明显减小；林芝站 (图 5c) 在500和250 hPa高度均方根误差略有所增大，其他高度误差都显著减小；昌都站 (图 5d) 除100 hPa高度，其他高度均方根误差均变小；巴塘站 (图 5e)、凉山站 (图 5f) 和宜宾站 (图 5i) 均方根误差变化程度不像其他站那么大，但误差也有所减小；红原站 (图 5g) 在中层误差减小十分显著；温江站 (图 5h) 在除400和150 hPa高度均方根误差变大外，其他高度均明显减小。总体而言，采用次网格地形重力波拖曳参数化方案之后，风场的均方根误差有较明显地变小，说明引入次网格地形重力波拖曳参数化方案对于这次个例风场模拟有正效果。

3.5.2 对降水的影响

以上分析表明次网格地形重力波拖曳对风场产生直接影响，而风场的变化必然会对降水造成一定的影响。为此，绘制了降水实况和两个试验模拟的降水分布 (图 6)；从图 6a实况可以看出，8日强降水中心主要位于四川省中北部成都、都江堰、彭州等地区，雨带呈西南—东北走向，24 h最大降水中心降水量超过150 mm。9日 (图 6b) 强降水中心仍位于成都、都江堰地区，24 h最大降水量也超过150 mm。对于0—24 h累计降水模拟，CTRL (图 6c) 试验模拟出两个强降水中心，一个降水中心位于四川西北部，降水极大值达到130 mm，另一个降水中心在四川省中部，最大降水量为110 mm，可见CTRL试验对于0—24 h累积降水的落区和强度的模拟均有较大偏差 (落区偏西、强度偏弱)，出现了较明显的漏报和错报现象；次网格地形重力波拖曳试验 (图 6e) 虽然仍模拟出两个强降水中心，但对于偏西的降水中心较CTRL试验降水量明显减小，而且，另一个一强降水中心位置和强度都更接近实况，强降水中心在成都、都江堰一带，最大降水中心降水量超过150 mm。对于24—48 h累积降水模拟而言，CTRL试验 (图 6d) 在四川地区完全没有模拟出强降水中心，有很明显的漏报现象；次网格地形重力波拖曳试验 (图 6f) 模拟出了降水中心的位置，但是降水强度较低，最大降水中心值为130 mm。

以上定性分析表明，引进次网格地形重力波拖曳参数化方案对降水模拟有一定改善，为了定量分析该参数化方案对降水模拟的影响，对这次个例次网格地形重力波拖曳和CTRL试验分别进行了0—24和24—48 h累积降水统计检验。检验结果如图 7所示，从CTRL和次网格地形重力波拖曳试验降水偏差可以看出，采用次网格地形重力波拖曳参数化方案后，0—24 h累积降水模拟 (图 7a) 除中雨模拟偏差略小外，其他等级降水都更接近于1；24—48 h累积降水模拟 (图 7b) 采用次网格地形重力波拖曳参数化方案后，从小雨至暴雨量级的降水预报偏差均有所减小。从0—24和24—48 h累积降水的TS评分 (图 7c、d) 看，除0—24 h中雨等级降水评分有所降低外，其他等级降水评分都有所提高或相当。

以上检验结果表明，总体来说，引入次网格地形重力波拖曳参数化方案之后改进了模式对这次个例降水的模拟效果。

4 批量试验统计检验结果

以上对风场和降水的检验都是针对一个个例，分析结论可能存在一定的偶然性，要评估引入次网格地形重力波拖曳参数化方案的总体性能，需要进行批量对比试验。为此，对2013年8月1—31日，进行了有、无次网格地形重力波拖曳方案的24 h模拟试验，模拟区域、分辨率以及物理过程方案的选取均与个例研究相同。

同样对前面提到的9个探空站进行了风场检验 (图 8)，由图可知，采用次网格地形重力波拖曳参数化方案后，拉萨站 (图 8a) 风场均方根误差在250 hPa变化不大，100 hPa误差变大约0.2 m/s, 其他高度误差明显减小；那曲站 (图 8b) 和林芝站 (图 8c) 除150 hPa均方根误差变大外，其他高度上误差均显著减小；昌都站 (图 8d) 在250和150 hPa高度上误差有明显减小；巴塘站 (图 8e)、凉山站 (图 8f) 和红原站 (图 8g) 均方根误差变化不太大；温江站 (图 8h) 在150和100 hPa高度上误差增大，其他高度上均有不同程度的减小；宜宾站 (图 8i) 在400—150 hPa误差减小明显，其他高度相差不大。总体来看，采用次网格地形重力波拖曳参数化方案后，各站风场的均方根误差都有明显减小，即风场模拟更接近实况。

由于青藏高原地形重力波的作用会影响到下游的天气，且青藏高原西侧和上空测站较为稀少，因此，对降水的检验主要考虑青藏高原东侧四川省附近。图 9为批量试验对降水的统计检验，有TS评分 (a)、漏报率 (b) 和空报率 (c)。从检验结果可知，采用次网格地形重力波拖曳参数化方案后，小雨等级降水的TS评分基本相当，其他等级降水的TS评分略有所提高，对于大雨等级降水评分提高比较明显。从漏报率看，除小雨等级降水略有增大外，其他等级降水都有一定程度的减小。从空报率来看，暴雨等级降水有明显减小，而大雨等级降水有明显的增大，小雨等级略有下降，中雨等级降水变化不大。综合来看，采用次网格地形重力波拖曳参数化方案对降水模拟具有正效果。

5 结论和讨论

针对目前对青藏高原大地形激发的重力波拖曳还不十分清楚的问题，在GRAPES_Meso模式中引入了KA95重力波方案和LM97阻塞拖曳方案，通过数值试验初步研究了青藏高原地区次网格地形重力波拖曳的一些基本参数，并对引入该方案后，数值模式的预报准确率进行了检验分析。得到的主要结论有：

(1) 沿30°N地形重力波拖曳的垂直分布显示，阻塞拖曳主要存在于模式层第1—5层 (第5层距地表高度约为417 m)，重力波拖曳主要存在于第5—10层 (第10层距地表高度约为1440 m)；从水平分布看，模式第3层以阻塞拖曳为主，主要位于青藏高原边缘地区，阻塞拖曳大值区沿喜马拉雅山脉走向和青藏高原东坡；模式第5层以重力波拖曳为主，主要位于青藏高原东部地区和云贵高原的北部边缘。

(2) 弗劳德数和气流绕流高度分析表明，在青藏高原喜马拉雅山脉和高原东部边缘地区，气流爬坡能力强，同时在这一地区绕流高度最高，绕流高度最高达第7层 (距地表高度约为757 m)；弗劳德数越大的地区绕流高度距离地表越高。

(3) 有、无采用次网格地形重力波拖曳参数化方案对比试验表明，采用该方案后，对于低层和高层地形重力波破碎有更准确的描述，地形重力波是向上垂直传播的。

(4) 个例和批量试验检验结果表明，数值模式采用次网格地形重力波拖曳参数化方案后，对于风场和降水模拟都有正效果，提高了模式预报的准确率。

本研究关于地形重力波拖曳的分析只是通过一个个例进行的，得到的分析结论也是初步的和探索性的。如果模式垂直分层或模式物理过程改变，研究结论可能也需要有相应的微调；但从批量数值试验的结果看，在数值模式中，加入地形重力波拖曳参数化方案对于提高模式的预报准确率是十分有益的。关于对青藏高原地形重力波拖曳的认识尚需进行更深入的研究和更多的观测事实来分析验证。