非静力地形重力波参数化方法技术

本发明专利技术公开了非静力地形重力波参数化方法，包含如下步骤：计算静力地形重力波地表动量通量和非静力效应对地形重力波动量通量的修正因子；计算非静力地形重力波地表动量通量；从位于地形高度之上的第一个模式层开始,进行如下步骤处理：计算地形重力波振幅；计算波动Richardson数以及地形重力波的饱和振幅和动量通量；若波动Richardson数大于临界值Ri

【技术实现步骤摘要】
非静力地形重力波参数化方法
本专利技术涉及数值模式中地形重力波参数化方法，尤其针对高分辨率数值模式中的非静力地形重力波参数化方法，属于大气科学研究领域。
技术介绍
地形重力波是数值模式的重要次网格物理过程，对于准确的数值天气预报和气候预测具有十分重要作用。传统数值模式水平分辨率较粗，因此次网格地形的水平尺度也较大，激发的重力波为静力平衡的波动(McFarlane1987)。近年来，随着数值模式的发展，模式的水平分辨率也不断提高，已经由百公里达到十公里量级，如ECMWF的IFS预报系统采用全球9km分辨率。在这种情况下，地形重力波的动力结构特征和动量传输将受到非静力效应的显著影响(Xueetal.2000)，因此为非静力地形重力波。对于非静力的地形重力波，通常无法求出解析解，因此无法得到非静力地形重力波的动量通量，给参数化带来很大挑战。为此，Smith(1980)提出采用二维快速傅里叶变换的方法对重力波进行求解，即数值求解。很显然，该数值方法无法用于地形重力波参数化方案。此后，Shutts(1998)等人采用物理空间的射线追踪方法(spatial-raytracing)来近似求解地形重力波。然而该方法仅能得到重力波的远场解，在地形正上方存在焦散点，波射线理论不再适应。为了克服这一难题，Broutmanetal.(2002)发展了谱空间中的波射线追踪方法，即Fourier-raytracing技术，有效避免了物理空间中地形上方的焦散点。然而该方法在浮力频率的转折点仍存在焦散点。近年来PulidoandRodas(2011)采用高阶的波射线近似方法(即高斯波束近似法)对地形重力波进行近似求解。与传统的波射线方法相比，该方法考虑了不同的波射线对中心波射线的贡献，因此能够有效避免传统射线追踪理论中的焦散点问题。利用该方法，Xuetal.(2017)研究了方向性切变对地形重力波动力结构和动量传输的影响，取得较好的效果。然而，高斯波束近似需要考虑不同的波射线贡献，计算量十分巨大，无法应用于业务化的地形重力波参数化。综上所述，为了满足高分辨率数值模式的发展需求，需要发展非静力地形重力波动量通量的计算方法。在此基础上，对传统的静力地形重力波参数化方案进行改进，发展非静力地形重力波参数化方案。
技术实现思路
本专利技术所要解决的技术问题是针对现有地形重力波参数化方案无法表征非静力效应，提出一种新的非静力地形重力波参数化方案，以满足高分辨率数值模式的发展需求，提高模式在复杂地形下的模拟预报能力。为了解决上述技术问题，本专利技术提出的非静力地形重力波参数化方法包含如下步骤：步骤1：计算静力地形重力波地表动量通量τs。步骤2：计算非静力效应对地形重力波动量通量的修正项α。步骤3：计算非静力地形重力波的地表动量通量τ0＝(1+α)τs。步骤4：从位于次网格地形高度之上的第一个模式层开始，进行如下步骤处理，步骤4.1：根据模式的水平风场Ui、Vi以及大气密度ρi和层结Ni,计算地形重力波振幅Hi。步骤4.2：根据地形重力波的振幅计算波动Richardson数Ri。步骤4.3：进行重力波破碎的判断。若波动Richardson数Ri小于临界值Ric，则地形重力波达到饱和，计算地形重力波的饱和振幅Hi_sat以及地形重力波在该模式层的动量通量τi；若波动Richardson数Ri大于临界值Ric，则地形重力波没有达到饱和，继续垂直传播，地形重力波在该模式层的动量通量与第i-1层相等，即τi＝τi-1。步骤4.4：重复步骤4.1至4.3，直至达到模式顶层。上述技术方案中，其中ρ0为地表空气密度，N0为地表大气层结，为次网格地形高度的平方，γ＝a/b为次网格地形各向异性，a、b分别为东西和南北方向的次网格地形宽度，U0、V0分别为东西和南北方向的地表水平风，χ为地表水平风方向，φ和为的方向和大小，其中k和l为东西和南北方向的水平波数。上述技术方案中，其中为水平Froude数。上述技术方案中，其中分别为重力波在地表和模式第i层的垂直波数，下标i表示模式的第i层。上述技术方案中，其中为平均气流的Richardson数，Uzi和Vzi分别为模式第i层东西风和南北风的垂直风切。上述技术方案中，上述技术方案中，上述技术方案中，Ric＝0.25。本专利技术的非静力地形重力波参数化方法具有以下几方面优点：1.能够表征非静力效应对地形重力波动量通量的影响。图1给出了各向同性地形激发的非静力重力波动量通量随着Fr数的变化。随着Fr数的增加，地形重力波动量通量逐渐减少。此外，由于Fr数与地形水平尺度有关，因此能够表示地形尺度变化所导致的地形重力波动量通量的变化，即尺度自适应。2.对于各向异性地形仍然具有较高的精度。图2给出了各向异性地形激发的非静力重力波动量通量与各向同性地形重力波动量通量的比值。随着Fr数的增加，各向异性对地形重力波动量通量的影响逐渐增加，但是误差不超过5％。因此，对于不同的各向异性地形，其地形重力波动量通量和各向同性的地形重力波动量通量十分类似。附图说明图1各向同性地形激发的非静力重力波动量通量随着Fr数的变化图。图2各向异性地形激发的非静力重力波动量通量与各向同性地形重力波动量通量的比值。图3全球地形(a)和次网格地形(b)的分布图。图4青藏高原西部地形重力波振幅,其中虚线为静力参数化方案，点线为非静力参数化方案和实线为饱和振幅的垂直分布。具体实施方式以2013年12月青藏高原西部地区的地形重力波为例，对本专利技术提出的非静力地形重力波参数化方案进行具体说明。第一步，地形和环境大气的预处理。根据欧洲中间尺度天气预报中心提供的2.5°再分析数据，可以得到包括青藏高原在内的全球地形和次网格地形分布，如图3所示。同时还可以得到2013年12月高原西部地区的风速以及浮力频率的垂直廓线。本步骤为人为在脱机的情况下处理，在模式系统中可直接跳过。第二步，计算静力地形重力波地表动量通量。该地形最大高度为5864米，因此计算该高度以下的平均风速、平均浮力频率以及平均密度，分别为U0＝6.8ms-1,V0＝2.1ms-1，N0＝0.0085s-1，ρ0＝1.12kgm-3。此外，由次网格地形数据可知，次网格地形东西和南北宽度分别为a＝20km,b＝16km，因此γ＝1.25。根据以上数据，最终计算得到静力地形重力波的地表动量通量为τs＝(0.000238,0.000817)Nm-1。第三步，计算非静力地形重力波动量通量。首先计算非静力因子α。根据数据可得水平Froude数为Fr＝0.311,代入α表达式可得α＝-0.011。因此非静地形重力波动量通量为τ0＝(1+α)τs＝(0.000212,0.000727)Nm-1。第四步，由于地形高度位于模式的第本文档来自技高网...

【技术保护点】
1.非静力地形重力波参数化方法，其特征在于：/n步骤1，计算静力地形重力波地表动量通量τ

【技术特征摘要】
1.非静力地形重力波参数化方法，其特征在于：
步骤1，计算静力地形重力波地表动量通量τs；
步骤2，计算非静力效应对地形重力波动量通量的修正项α；
步骤3，计算非静力地形重力波的地表动量通量τ0＝(1+α)τs；
步骤4，从位于次网格地形高度之上的第一个模式层开始，进行如下步骤处理；
步骤4.1，根据模式的水平风场Ui、Vi以及大气密度ρi和层结Ni，计算地形重力波振幅Hi；
步骤4.2，根据地形重力波的振幅计算波动Richardson数Ri；
步骤4.3，若波动Richardson数Ri小于临界值Ric，计算地形重力波的饱和振幅Hi_sat以及地形重力波在该模式层的动量通量τi；若波动Richardson数Ri大于临界值Ric，继续垂直传播，地形重力波在该模式层的动量通量与第i-1层相等，即τi＝τi-1；
步骤4.4：重复步骤4.1至4.3，直至达到模式顶层。


2.如权利要求1所述的非静力地形重力波参数化方法，其特征在于：



其中ρ0为地表空气密度，N0为地表大气层结，为次网格地形高度的平方，γ＝a/b为次网格地形各向异性，a、b分别为东西和南北方向的次网格地形宽度，U0、V0...

【专利技术属性】
技术研发人员：徐昕，李闰秋，张蓉蓉，
申请(专利权)人：南京大学，
类型：发明
国别省市：江苏;32