引言

自20世纪70年代中期以来, 已有不少研究数值模式静力与非静力差别的数值试验和理论^[1-2]。Daley^[3]研究发现, 用静力与非静力数值模式计算球面大气波的正规模, 对于纬向波数大于400的大气波动, 两者相去甚远, 从而表明不宜用静力模式描述小于100 km的大气波动。Dudhia^[4]设计建立非静力全可压、具有较为完全物理过程的中尺度数值模式MM5;类似作为业务化的中尺度或全球/区域通用非静力数值模式还有:ALADIN^[5], ARPS^[6], JMA_MRI^[7], UKMO^[8], RAMS^[9], 中国的GRAPES_ meso^[10]; 及近年由美国多个研究单位和大学联合开发的新一代细网格(水平分辨率达到1~10 km)中尺度数值模式WRF (Weather Research and Forecasting)。

现阶段非静力模式均须引入所谓的"静压"参考大气(参考大气仅为高度函数, 也称为"参考大气廓线"), 同时对气压、气温场作"扰动"线性化处理, 以便扣除垂直(运动方程中)气压梯度力与重力有着静力平衡约束关系的参考大气重量, 使之成为所谓"非静压扰动平衡, 从而达到提高计算垂直气压梯度力与上游点垂直位移精度的目的。但参考大气一般只能取为初始大气平均"、"Lorenz参考大气"^[11]、"等温大气"或气候平均"等, 这些"参考大气廓线"均与实际(模式)三维大气内在的、变化的静力平衡态有差距。廖洞贤^[12-13]指出, 现在的非静力模式所计算的非静压扰动一般约占整个气压场10%或稍小, 而模式大气的非静压扰动与"静压"参考大气之比应约为1:99, 则引入"参考大气廓线"之后仍然存在较大的误差, 从而否定模式大气中存在"定常的"参考大气廓线, 进而指出, 静压参考大气只能是"时变的"时空函数。

研究数值模式非静力全可压动力框架, 有一个国内外认可、有效、可对比的理想试验方法, 即密度流试验, 用以检验数值模式对于高度非线性流的模拟水平。1990年9月在美国大型计算中心(National Center for Super- computing Applications, NCSA)举行"非线性流计算方法研讨会", 提出了密度流试验的benchmark参考解^[14], 认为各种数值计算方法的结果应该收敛于benchmark参考解, 包括它的基本物理特征及大部分细节, 则其成为检验新格式非静力模式数值解的一致性、收敛性及稳定性的参考解, 但奇怪的是, 当时的高阶格式方案和谱方法做出的最好结果都是无效、不可信的, 比较而言, "单调格式"(指以一阶中央差为代表的各种欧拉线性格式)却可以给出较为一致、连续的结果, 但是, 单调格式即便采用高时空分辨率, 也总是存在对波动位相和振幅描述显著的误差。我国研究开发的GRAPES全球/区域统一数值模式, 设计为非静力全可压动力框架和采用半隐式-半拉格朗日积分方案和Lagrangian三次函数插值(双线性插值), 垂直方向是采用Charney-Philips跳点格式, 由杨学胜等^[15]完成了密度流试验。

本文将采用高阶(二阶精度)、非线性(二阶可导)样条格式和显式准拉格朗日积分方案做密度流试验, 并认为, 若是采用适当的平滑方案(未来应可采用三次样条光顺泛函数, 借以消除局部"多余拐点"和不连续"尖点"或"环绕"), 则准拉格朗日法样条格式完全可以优于单调格式, 因在数学上, "样条格式+局部平滑"是对波动(位相和振幅)作"二阶可导(三次连续)"、"收敛"和最优"描述^[16](中央差斜率和曲率只是三次样条斜率和曲率作系数为1/3的"三点平滑")。众所周知, 样条格式既为差分又兼积分格式, 则通过做气压、气温场双三次曲面拟合就能计算(已隐含)全球大气质量和球面大气能量, 且做静力方程三次样条并积分, 就能计算"时空函数"静压参考大气。故我们认为, 最终应能实现全球/区域(多重嵌套)统一样条格式、静力/非静力(全可压)数值模式动力框架。下面列出当前研究非静力数值模式的几个基本格局(可与本文密度流试验相比较):

(1) Charney-Phillips垂直跳点格式难以设计成为二阶差分精度;

(2) 二时间层迭代求上游点路径, 考虑两个时间层之间有风速变化(有加速度), 是用平均风速求路径, 本文考虑按一时间层(当前时刻)风速和加速度迭代求上游点路径, 二者均符合牛顿运动定律(后者完全二阶时空精度高于前者);

(3) GRAPES模式对Exner指数气压和位温做预报, 本文考虑对(自然对数)气压、气温/位温做预报, 但在对气温/位温场、尤其对垂直气温/位温场作平滑时, 需用到位温守恒原理;

(4) 引入参考大气, 用以提高计算垂直气压梯度力(加速度)及位移精度, 本文考虑计算出时变的静压参考大气, 准确(二阶空间精度)计算垂直加速度, 更无须引入参考大气廓线;

(5) 在引入参考大气之后, 又将风变率(单位质量空气广义牛顿力或加速度)分为线性项和非线性项两个部分, 做非时间中央差积分, 其中非线性项采用二时间层外推与三时间层平滑, 本文考虑将风变率(广义牛顿力=气压梯度力+地转偏向力+曲率力+摩擦力+重力)当作一个不可分割的整体, 通过三次样条函数也拟合成为一个二阶可导变量场, 且做一时间层"匀加速"时间积分;

(6) 采用SI-SL积分方案, 必须求解庞大计算量的(三维气压扰动)广义共轭余差Helmholtz方程, 本文研究样条格式显式-准拉格朗日积分方案, 求得二阶时空精度上游点压、温、湿、风及加速度与三维(位移)散度场, 其计算量较小, 且三次样条函数适合做并行计算。

1 准拉格朗日非静力全可压动力框架 1.1 原始大气运动方程

为做理想初值场密度流试验, 采用如下的局地直角坐标系原始大气运动方程^[17](不考虑地形, 仅列出干绝热、干过程压、温、湿、风预报方程):

$ \frac{{{\text{d}}u}}{{{\text{d}}t}} = - RT\frac{{\delta \ln p}}{{\delta x}} - fv - \tilde fw + {F_u} = {a_u} $

(1)

$ \frac{{{\text{d}}v}}{{{\text{d}}t}} = - RT\frac{{\delta \ln p}}{{\delta y}} - fu + {F_v} = {a_v} $

(2)

$ \frac{{{\text{d}}w}}{{{\text{d}}t}} = - RT\frac{{\delta \ln p}}{{\delta z}} - g - \tilde fu + {F_w} = {a_w} $

(3)

$ \frac{{{\text{d}}\ln p}}{{{\text{d}}t}} = \frac{{ - 1}}{{1 - k}}\nabla \cdot V = {a_p} $

(4)

$ \frac{{{\text{d}}\ln T}}{{{\text{d}}t}} = \frac{{ - k}}{{1 - k}}\nabla \cdot V = {a_T} $

(5)

$ q = \frac{{{\text{d}}q}}{{{\text{d}}t}} \equiv 0 $

(6)

1.2 样条格式二阶时空离散预报方程

可按文献[18]准拉格朗日法二阶时空离散显式预报方程预报(u，v，w)风场, 而按隐式预报方程预报气压、气温场, 则相应(1)-(6)式给出如下密度流试验的非静力全可压动力框架(定义Θ=lnp):

$ {u^{t + \Delta t}} = {{\dddot u}^t} + \dddot a_u^t\Delta t $

(7)

$ {v^{t + \Delta t}} = {{\dddot v}^t} + \dddot a_v^t\Delta t $

(8)

$ {w^{t + \Delta t}} = {{\dddot w}^t} + \dddot a_w^t\Delta t $

(9)

$ \begin{gathered} {p^{t + \Delta t}} = {{\dddot p}^t}\exp \left( { - \frac{{\Delta t}}{{1 - k}}\Delta \cdot \bar V} \right)或 \hfill \\ {\Theta ^{t + \Delta t}} = {{\dddot \Theta }^t} - \frac{{\Delta t}}{{1 - k}}\Delta \cdot \bar V \hfill \\ \end{gathered} $

(10)

$ {T^{t + \Delta t}} = {{\dddot T}^t}\exp \left( { - \frac{{k\Delta t}}{{1 - k}}\Delta \cdot \bar V} \right) $

(11)

$ {q^{t + \Delta t}} \equiv 0 $

(12)

上面(10)式表明, 原始大气运动方程既允许对气压(p), 又允许对自然对数气压(ln p)做预报(本文的密度流试验就是对ln p做预报)。

上面气压、气温预报方程(10)、(11)式均含一个时间步长(Δt)离散、平均三维散度$\Delta \cdot \bar V$。对$\Delta \cdot \bar V$是先取其平均风速$\bar V = \left({\bar u, \bar v, \bar w} \right) = \left({\frac{{\Delta x}}{{\Delta t}}, \frac{{\Delta y}}{{\Delta t}}, \frac{{\Delta z}}{{\Delta t}}} \right)$, $\left({\Delta x, \Delta y, \Delta z} \right)$为三维位移, 再对$\left({\Delta x, \Delta y, \Delta z} \right)$分别做(x, y, z)方向三次样条拟合, 则容易求得(10)式和(11)式(自然对数)气压和气温"样条格式"隐式变率:$\Delta t\nabla \cdot \bar V = \Delta {x^x} + \Delta {y^y} + \Delta {z^z}$, 即是用它代替隐式预报方程中上游点"起点变率+终点变率"平均值:

$ \frac{{a_p^{t + \Delta t} + \dddot a_p^t}}{2} = \frac{{ - \nabla \cdot \bar V}}{{1 - k}} $

(13)

$ \frac{{a_T^{t + \Delta t} + \dddot a_T^t}}{2} = \frac{{ - k\nabla \cdot \bar V}}{{1 - k}} $

(14)

上述预报方程通过"显式"迭代插值求三维上游点, 既是对气块作样条格式三维位移"(三次运动"质量平流)描述, 同时, 又是"隐式"求其三维(位移)散度场。这里物理机制是, 上游点气块受到三维散度场作用(所谓"全可压"), 其压、温(密度)发生绝热变化, (此时还应考虑气块水汽是否达到饱和与是否发生凝结降水)。故上述预报方程, 包含了从声波(声波为压缩波, 声波波速约为330 m·s^-1)到重力波(重力内波波速约为30 m·s^-1)等快波。所以, 上述动力框架与密度流试验, 只能将就声波波速而采用极高的时空分辨率。

2 两种方法计算垂直气压梯度力

准确计算垂直气压梯度力, 关系到模式"非静力全可压"动力框架精确性。

传统方法是引入一个定常的(干)参考大气廓线"(参考廓线法"), 其参考气压$\tilde p\left( z \right)$和参考气温$\tilde T\left( z \right)$应满足:

$ \frac{{\partial \ln \tilde p}}{{\partial z}} = \frac{g}{{{R_d}\tilde T}}\left( {静力平衡方程} \right) $

(15)

本文方法是考虑气柱隐含一个"时空函数"静压场，${\bar p}$(x，y，z，t) "(分离静压法"), 其满足:

$ \frac{{\partial \ln \bar p}}{{\partial z}} = \frac{g}{{RT}}\left( {{\text{静力平衡方程}}} \right) $

(16)

或更精确地满足球面坐标含压、温、湿、风的"静力平衡方程"^[18]:

$ \frac{{\partial \ln \bar p}}{{\partial z}} = \frac{{g - \tilde fu - \left( {{u^2} + {v^2}} \right)/{r_0}}}{{RT}} $

(17)

上式r₀为地球平均半径, 实际应用时, 地球半径与重力加速度(g)可以不为常数。

对(16)式右边(含温、湿及重力场), 或对(17)式右边(含温、湿、风及地球半径与重力场), 作气柱的三次样条拟合(顶/底层取前/后差分边界), 并设顶层(K层) $\bar p = {p_k}$为常数层, 其物理意义是, 整个模式顶层之上大气质量守恒, 则可以从模式顶层向底层积分, 求得各层的${\bar p}$, ${\bar p}$是气柱"时变的"分层重量分布。同理, 对(15)式右边作气柱三次样条拟合, 求得各层的$\tilde p$, ${\bar p}$是"参考大气廓线"分层重量分布。

对应于(15)式$\tilde p$和(16)式(或(17)式) ${\bar p}$, 分别设$p = \tilde p \cdot \tilde p'$和$p = \bar p \cdot \bar p'$, 容易分别求得它们各层${\tilde p'}$和${\bar p'}$, 再对各层$\ln \tilde p'$和$\ln \bar p'$作垂直三次样条, 求得${\left({\ln \tilde p'} \right)^z}$和${\left({\ln \bar p'} \right)^z}$, 但顶、底层均取为静力平衡边界, 以使得在顶、底层(边界)上${\left({\frac{{\delta \ln \tilde p'}}{{\delta z}}} \right)_k} = {\left({\frac{{\delta \ln \tilde p'}}{{\delta z}}} \right)_0} \equiv 0$和${\left({\frac{{\delta \ln \bar p'}}{{\delta z}}} \right)_k} = {\left({\frac{{\delta \ln \bar p'}}{{\delta z}}} \right)_0} \equiv 0$, 即使得它们的垂直加速度为(a_w)_K≡0和(a_w)₀ ≡0, 且密度流试验将设定摩擦作用F_w≡0, 则求得它们(球面坐标w方程)^[18]垂直加速度a_w (w变率)分别为:

$ \begin{gathered} {a_w} = - RT\frac{{\delta \ln \tilde p'}}{{\delta z}} - g\left( {1 - \frac{{RT}}{{{R_d}\tilde T}}} \right) + \tilde fu + \frac{{{u^2} + {v^2}}}{{{r_0}}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\; \approx - RT{\left( {\ln \tilde p'} \right)^z} - g\left( {1 - \frac{T}{{\tilde T}}} \right) \hfill \\ \end{gathered} $

(18)

$ {a_w} = - RT\frac{{\delta \ln \bar p'}}{{\delta z}}\; \approx - RT{\left( {\ln \bar p'} \right)^z} $

(19)

上面(18)式(略去了曲率项和离心力项2个小项)和(19)式分别为参考廓线法和分离静压法样条格式求垂直加速度a_w, 本文密度流试验采用(18)式计算a_w时, 因采用局地坐标而舍去了曲率项和离心力项, 采用(19)式计算a_w时, 实际也采用局地坐标(16)式而不是球面坐标(17)式。

上面(19)式表明, 分离静压法已从模式大气非静力气压场完全分离出时变的静压场, 从而使得计算垂直气压梯度力无须再引入模式大气中并不存在的"参考大气廓线", 这将能够大大简化预报方程, 无须将气压、气温变量再人为分为"基本量+扰动量", 同时能够提高计算垂直气压梯度力与位移的精度, 这里扣除的静压场, 甚至能够通过给定"时空函数"重力场, 而考虑到地球椭球性和月球引潮力。

还应该指出, 很容易找到适合做密度流试验的参考大气廓线", 即是密度流试验未受扰动的初值场(图 1-2), 实际强对流过程多种多样, 很难找到适合预报对象的"参考大气廓线", 且不得不将气压、气温变量作所谓"线性化"处理, 这样做既违背牛顿运动定律, 也使得预报方程复杂化。

3 密度流试验

密度流试验是在局地直角坐标系与均匀网格上完成。密度流试验将检验采用样条格式计算"非静力全可压"气块路径的可行性, 即要检验压、温、风"三次运动"平流、弯流(及扭流)的基本特性与精确性。

密度流试验是一个二维(x, z)理想场试验, 初值场仍取为三维模式大气, 但仅对(y)方向的中间垂直截面预报点做时间积分, 并在初值场及每一时间步, 均对(y)方向其它预报点赋予相同的预报变量值(这时P^y=P^xy≡0, 未考虑(y)方向平流、弯流和扭流)。典型密度流试验分辨率取为Δx=Δy=Δz=100 m, 无地形, (x, y, z)方向分别取为(0:512)、(-4:4)和(0:53)格点分布, (x, y)方向均取周期三次样条边界, 而z方向的顶、底层, 气压场, 包括扰动量lnp′场, 均可取为静力平衡边界, 这样使得顶、底层w=a w≡0(刚体), 但气温和水平风场等, 在顶、底层均取为前、后差分边界, 密度流试验的benchmark参考解要求时间步长取为0.1 s, 积分900 s, 而本文采用"分离静压法"做密度流试验时, 还增加将时间步长取为0.05 s, 且增加延长积分到1 200 s。

未扰动初值场为干静止大气, 即q=u=v=w=0, 并地面层气压为1 000 hPa, 且整个模式大气初值位温(θ)场均为300 K, 则地面层气温为300 K。

但于未扰动初值场的中心放置一个圆形(截面)的冷涌(图 1), 即设:

$ \begin{gathered} \Delta \theta \frac{{ - 15}}{2}\left[ {\cos \left( {\pi \cdot L} \right) + 1} \right], \hfill \\ L = \sqrt {{{\left( {\frac{{{x_i} - {x_c}}}{{{r_x}}}} \right)}^2} + {{\left( {\frac{{{z_k} - {z_c}}}{{{r_z}}}} \right)}^2}} \leqslant 1, \hfill \\ i = 0, 1, \cdot \cdot \cdot , 512;k = 0, 1, \cdot \cdot \cdot , 53 \hfill \\ \end{gathered} $

(20)

这里, 取(x_c, z_c)=(256×100 m, 30×100 m)为冷涌中心点, (r_x, r_z)=(40×100 m, 20×100 m)为冷涌(x, z)方向半径。则在冷涌中心点:L=0, Δθ=-15K。并按静力平衡方程和位温守恒Poisson方程, 不难求得扰动初值场各层的气压和气温分布(图 1), 因冷涌植入, 其正下方的地面气压也有所增加, 达到1 013.21 hPa。

故密度流试验的预报变量为自然对数气压lnp、气温T (位温θ仅是诊断物理量)、风场u和w。

参考廓线法:实际取${a_w} = - {R_d}T\frac{{\partial \ln \tilde p'}}{{\partial z}} - g\left({1 - \frac{T}{{\tilde T}}} \right)$, 因各层$\tilde p' = p/\tilde p$(或$\ln \tilde p' = \ln p - \ln \tilde p$), 参考大气廓线$\tilde p\left(z \right)$和$\tilde T\left(z \right)$取为未扰动初值气压、气温场(均保持为常数)。例如, 扰动初值场各层$p = \tilde p$, 则各层$\tilde p' = 1$(或$\ln \tilde p' = 0$), 则有${a_w} = - g\left({1 - \frac{T}{{\tilde T}}} \right)$(图 2)。

分离静压法:实际取${a_w} = - {R_d}T\frac{{\partial \ln \bar p'}}{{\partial z}}$, 各层$\bar p' = p/\bar p$(或$\ln \bar p' = \ln p - \ln \bar p$), ${\bar p}$(x, y, z, t)为时变(时空函数)气压场。例如, 扰动初值场各层$p = \bar p$, 则有${a_w} = - {R_d}T\frac{{\partial \ln \left({\tilde p/\bar p} \right)}}{{\partial z}} = - {R_d}T\left({\frac{{\delta \ln \tilde p}}{{\delta z}} - \frac{{\delta \ln \bar p}}{{\delta z}}} \right) = - g\left({1 - \frac{T}{{\tilde T}}} \right)$(将(15)和(16)式代入), 这里证明, 两种方法的初值广义加速度相等(图 2)。

高度非线性的密度流试验, 揭示整个"冷涌下沉→底层冷空气堆积→冷锋前沿Kelvin-Helmholtz水平风切变形成→不稳定涡旋形成"演变过程, 且为"声波+重力波"传播过程。

分离静压法密度流试验, 其干绝热过程时间积分演变(图 3a-c、图 4a-c、图 5和图 6):在垂直气压梯度力(负浮力)作用下, 冷气流逐渐加速下沉, 形成下沉辐散气流, 在冷涌到达底层之前, 冷涌正下方地面气压从1 013 hPa迅速下降, 一度降到约1 002 hPa, 与地面气压降低相对应, 冷涌正下方近地面层900 m始终为下沉运动, 垂直风速在积分200 s时达到约-14 m·s^-1, 在冲击到达底层之后冷涌开始堆积。冷涌冲击过程分为正向压迫与反向回弹(所谓"全可压"或完全弹性"), 首次冲击到底层时, 冷涌正下方地面气压一度又增加至1 013 hPa, 首次反弹时, 地面气压又回到约1 005 hPa, 冲击波振幅约为7~8 hPa, 冲击波过程约在30 s内完成, 之后进入多个振幅约为3 hPa或更小振幅的"完全弹性"小波动, 小波动振幅约为2~3 hPa, 可见冲击波和小波动均为重力波, 下一次冲击与上一次冲击相间约150 s, 几经循环, 且冲击波与小波动都向着逐渐减压(向地面平均气压1 000 hPa)方向发展。

首次冷涌冲击波及冷空气堆积在地面形成了较强的水平气压梯度(力), 同时, 冷空气开始"一分为二"、两两相背在水平方向运动, 并形成"冷锋"(图 3- 4仅显示冷锋向x方向对称运动的正方向右侧一支)。此外, 与冷涌正下方地面气压变化相似, 在冷涌右左两侧10 km的地面气压也呈冲击波与小波动演变(图 5), 但后者冲击波振幅约为3~4 hPa, 小波动振幅约为1 hPa, 两次冲击波相间约75 s, 与前者比较, 后者振幅较小, 但频率较快。同时(图 7), 冷涌正下方地面水平风仅表现为声波震动, 其振幅约为±0.002 m × s^-1, 而在冷涌右侧/左侧10 km地面水平风表现为重力波冷锋运动特征, 在冷锋过境前水平风速逐渐增大, 在积分600 s冷锋过境时最大水平风速达到约23 m × s^-1/-23 m × s^-1, 不难看出, 整个水平风演变都包含了重力波与声波传播, 其中周期性叠加了最大振幅约为±1 m × s^-1的重力波, 而与地面水平风相对应, 冷涌右侧/左侧10 km的近地面层900 m垂直风(图 6)在冷锋经过时最大上升运动也达到约7.5 m × s^-1, 此后上升风速迅速减小, 一度还变为约-1 m × s^-1的下沉运动, 但在副冷锋经过时又再度变为上升运动且达到约5.5 m × s^-1。

分离静压法密度流试验结果表明(图 3a-c和图 4a-c):时间积分300 s之后, 冷涌主体到达底层, 并且在冷锋前方形成较强的水平风垂直切变, 达到Kel- vin-Helmholtz切变不稳定条件, 则第1个涡旋形成(图 3a、图 4a); 积分600 s, 第1个涡旋迅速发展增强, 并在减弱冷锋前方形成第2个涡旋(图 3b、图 4 a); 至积分900 s, 第1个涡旋发展成一个圆涡, 第2个涡旋也进一步发展(图 3c和4c), 但与Straka等^[14]给出单调格式相同分辨率的benchmark参考解相比较, 弱冷锋前方第3个涡旋尚未能生成, 但之后(积分1 200 s的图 8和图 9)还是完整地形成了(两者存在差别, 除它们采用了不同格式主要原因之外, 还应是它们采用了不同平滑方案和不同侧边界所造成)。

参考廓线法密度流试验(积分900 s)结果(图 3d、图 4d)与上述分离静压法试验(积分900 s)结果(图 3c、图 4c)两相比较, 它们模拟的气压和气温(位温)场均非常相似, 仅后者模拟的冷锋移速稍快于前者。

以上两种方法密度流试验结果表明, 本文分离静压法:无须引入"静压"参考大气廓线, 仅通过做气柱静力方程三次样条拟合与垂直积分, 就能够分离出时变的静压场, 且能够准确计算垂直气压梯度力。

4 密度流试验比较与分析

因上述两种计算垂直气压梯度力方法的密度流试验结果相近、且都接近于benchmark参考解, 从而表明, 本文非静力全可压动力框架"样条格式二阶时空余差预报方程+显式-准拉格朗日积分方案", 能够模拟高度非线性密度流的基本特征, 但若仔细比较其与benchmark参考解的差异^{[14, 19-21]}, 可作如下分析(仅对分离静压法密度流试验与结果作叙述):

4.1 样条格式与中央差关系

非线性样条格式是以"3点"为主但关于全局(全部格点)的"二阶可导"三次格式, 即需已知4个(以上)格点值和2个端点边界条件才可以做三次样条拟合^[16]。

对于离散等距ΔX变量场P_j, j=0, ± 1, ± 2, ..., ±J, 三次样条拟合斜率m_j和曲率M_j与一阶中央差${{\bar m}_j}$和二阶中央差${{\bar M}_j}$有如下数学关系(样条格式与中央差数学关系)^[16]:

$ {{\bar m}_j} = {m_j} + \frac{{{m_{j + 1}} - 2{m_j} + {m_{j - 1}}}}{6} = \frac{{{P_{j + 1}} - {P_{j - 1}}}}{{2\Delta x}} $

(21)

$ {{\bar M}_j} = {M_j} + \frac{{{M_{j + 1}} - 2{M_j} + {M_{j - 1}}}}{6} = \frac{{{P_{j + 1}} - 2{P_j} + {P_{j - 1}}}}{{\Delta {x^2}}} $

(22)

两式表明, 中央差斜率和曲率, 分别是三次样条斜率和曲率作系数为1/3的"三点平滑"。但样条格式曲率M_j的"线性部分"就是二阶中央差${{\bar M}_j}$^[16], 而样条格式斜率m_j由曲率M_j决定(二者不独立)。

4.2 样条格式边界

三次样条存在多种数学边界, 但需给定①边界斜率(一阶导数)或②边界曲率(二阶导数), 且存在周期三次样条边界(周期三次样条无边界)。

本文密度流试验顶/底层气压场取为静力平衡边界, 其它如气温和风场等均取为前/后差分边界, 而(x)方向侧边界均取为周期三次样条边界。由于密度流为一个(x)方向对称运动, 对于p、T、w场, 在两个侧边界上, 总有P_-J≡P_J和P_-J+1≡P_J-1, 这时有:

$ {{\bar M}_{ - J}} = \frac{{{P_{ - J + 1}} - 2{P_{ - J}} + {P_{J - 1}}}}{{\Delta {x^2}}} = \frac{{{P_{ - J + 1}} - 2{P_J} + {P_{J - 1}}}}{{\Delta {x^2}}} = {{\bar M}_j} $

(23)

而对于u场, 在两个侧边界上, 总有u_-J≡-u_J和u_-J+1 ≡-u_J-1, 这时有:

$ {{\bar M}_1} = \frac{{{u_{ - J + 1}} - 2{u_{ - J}} + {u_{J - 1}}}}{{\Delta {x^2}}} = \frac{{{u_{ - J + 1}} - 2{u_J} + {u_{J - 1}}}}{{\Delta {x^2}}} = - {{\bar M}_j} $

(24)

以上表明, 在密度流试验两个侧边界上, 它们p、T、w场曲率总相等, 而它们斜率(样条格式斜率由曲率确定)总"大小相等、方向相反", 但它们u场曲率总"大小相等、方向相反", 而它们斜率总相等, 即是它们在作关于(x)方向对称的相背或相向运动(在周期边界上为作相向的"对撞"运动)。

本文密度流试验采用样条格式垂直边界:顶、底层气压场为静力平衡边界(在顶、底层边界上广义加速度a_w≡0), 其它变量场(暂)为前/后差分边界, 而侧边界均取为周期三次样条边界, 将来还可以(试验)采用其它多种样条格式垂直边界或侧边界。

4.3 样条格式斜率、曲率和挠率

与单调格式相比较, 样条格式可以考虑变量场斜率、曲率和挠率拟合(二维密度流试验忽略了水平挠率), 其为完全二阶空间精度。因同元变量可有不同的物理含意与数学表述及空间形状, 例如(图 1), 初值扰动位温场呈椭圆形, 而相应的扰动气温场呈倒"Ω"形, 它们具有不一样的斜率和曲率。本文通过将密度流试验预报变量由(lnp, T, u, w)改变为(lnp, θ, u, w)时, 即是将预报气温变为预报位温时, 两个分离静压法密度流试验(积分1 200 s)结果(图 8a、b、图 9a、b)便有所不同。同理, 若是将预报变量由自然对数气压(lnp)改变为气压(p), 也将得到不完全一样的结果(暂未进行对比试验)。

三维散度场是预报气压、气温(及水汽相变)场的变率, 在先求得样条格式二阶时空精度的上游点三维运动路径之后, 再采用样条格式求得三维位移隐含的、一个时间步长平均三维散度场, 则本文密度流试验是采用了"变加速"隐式预报方程做气压、气温场预报。

4.4 样条格式与平滑

以解析的"二阶可导"样条格式模拟高度非线性的密度流试验, 其拟合气压、气温场、尤其拟合风场, 也会出现不稳定:即出现三次样条拟合之尖点或环绕^[16], 必需要有适当的滤波过程。本文尚未给予滤波, 仅对模拟过程做了简单的平滑处理, 若不做任何平滑, 则在模拟过程之中, 将出现计算不稳定, 造成模拟失败。

在密度流试验过程中, 每三个时间步(0.3 s)做一次垂直风(w)场在垂直方向、系数为1/3的三点平滑, 同时做一次气压(lnp)场在纬向、系数为1/2的三点平滑, 且在平滑气压场时, 对气温(T)场做相应的Poisson方程"绝热变温"; 又每900时间步(90 s)做一次气温场分别在纬向和垂直方向、系数均为1/2的"位温保守三点平滑。

本文密度流试验所采用的"位温保守"气温场纬向(i=0, ±1, ±2, ..., ±256)与垂直方向(k=0, 1, 2, ..., 53)三点平滑(${{\bar T}_i}$、${{\bar T}_k}$, c为平滑系数)分别为:

$ {{\bar T}_i} = {T_i} + \frac{c}{2}\left( {{T_{i + 1}}\frac{{p_i^k}}{{p_{i + 1}^k}} + {T_{i - 1}}\frac{{p_i^k}}{{p_{i - 1}^k}} - 2{T_i}} \right) $

(25)

$ {{\bar T}_k} = {T_k} + \frac{c}{2}\left( {{T_{k + 1}}\frac{{p_k^k}}{{p_{k + 1}^k}} + {T_{k - 1}}\frac{{p_k^k}}{{p_{k - 1}^k}} - 2{T_k}} \right) $

(26)

上述(25)和(26)式也可以用作预报变量位温(θ)场平滑, 但同时应该做相应的气温(T)场"诊断"平滑。

本文密度流试验结果已经表明, 做"位温保守"气温场平滑是一种可取的平滑方法。将来还可考虑按变量场三次样条曲率去判断大气运动必然存在的临界不稳定点(域), 以便做变量场单点或局域平滑, 并且在平滑之后, 容易重新求得平滑点(域)二阶可导"补丁"。而将来研究应用三次样条光顺泛函数, 借以消除局部的多余拐点、不连续尖点或环绕, 尤显重要。(采用不同边界与平滑方案是造成本文密度流试验结果与benchmark参考解"细节"有差别的另两个重要原因)。

4.5 样条格式解析性

样条格式具备二阶可导"解析性", 样条格式将压、温、(湿)场、以及旋转地球上的风场与单位质量空气所受广义牛顿力(加速度(a_u, a_v, a_w))场都拟合成为二阶可导, 使得准拉格朗日法计算出解析的上游点与伴随上游点的广义加速度, 即样条格式可以将(a_u, a_v, a_w)拟合为二阶可导的(当作不可分割的)连续场, 从而能(且仅能)采用"匀加速"显式预报方程做风场预报, (而隐式和半隐式积分方案需要将广义加速度中的曲率项和含垂直速度的非线性小项舍去, 从而便于隐式求解水平风场, 但舍去多个小项, 将使得广义加速度失于连续性, 且有悖于牛顿运动定律)。

4.6 样条格式收敛性

样条格式具备二阶可导"收敛性", 其数学表述是:"三次样条及其一阶、二阶导数收敛于原函数及其一阶、二阶导数", 文献[18]中Rossby-Haurwitz波试验, 采用不同时间步长, 试验结果已经证明样条格式具有收敛性。本文也将密度流试验时间步长由0.1 s改变为0.05 s, 两个(参考廓线法)试验做相同时间间隔的平滑(平滑总次数相同), 它们(积分1 200 s)的试验结果(图 8a、c、图 9a、c)却较为相近, 其收敛性似"不明显", 我们认为, 这或许是全局平滑"歪曲"作用远大于(掩盖了)其收敛性。

4.7 样条格式最优性

样条格式具备二阶可导"最优性", 其数学表述是:"三次样条二阶导数对原函数二阶导数为最佳逼近", 而三次样条一阶导数(斜率)由二阶导数(曲率)确定, 综合上述密度流试验结果及比较分析, 我们认为, 应用样条格式做非静力全可压数值模式动力框架, 最终应该优于(不会输于)单调格式。

5 结论和讨论

(1) 研究建立样条格式非静力全可压动力框架, 因三次样条函数是二阶空间精度"差分+积分"格式, 其以"二阶可导"拟合压、温、湿、风及广义牛顿力场, 从而能实现样条格式显式-准拉格朗日积分方案:采用"匀加速"显式预报方程做风场预报, 采用"变加速隐式预报方程做压、温场预报, 因压、温场变率是三维(位移)散度场, 故整个"显式+隐式"积分方案, 本质上只是显式方案。

(2) 样条格式密度流试验:采用三次样条函数(线、面、体)去拟合各个变量场, 拟合了变量场斜率(平流)、曲率(弯流)和挠率(扭流), 但密度流试验为二维试验, 未用到变量场水平挠率, 但将来发展非静力全可压数值模式动力框架, 仍需要对三维变量场做三次曲体拟合, 即需发展所谓"曲体元"模式, 这对于预报强对流天气系统(如下击暴流)是必要的。

(3) 本文的高度非线性密度流试验结果, 初步证明样条格式显式-准拉格朗日积分方案动力框架具备模拟精细尺度非线性流及其瞬变特征的能力, 验证样条格式用做格点数值模式动力框架的(与其它单调格式)一致性、(二阶精度)精确性及程序正确性。

(4) 密度流试验表明, 因三次样条函数是一种(二阶空间精度)积分格式, 通过做静力平衡方程三次样条拟合, 可以从非静力气压场分离出时变静压场(时空函数静压场), 从而无须引入参考大气, 便能准确计算垂直气压梯度力与垂直位移, 而传统方法是给定一外在的参考大气廓线, 一般只由"等温大气"或"平均大气代替, 同时, 对应地将气压、气温场人为分为"基本量+扰动量"两个部分, 目的也为提高计算垂直气压梯度力的精度, 而在非静力模式中, 实际很难找到合适的(像密度流试验"未经扰动"初始场一样的)参考大气廓线, 可以说, 对于模式大气, 引入参考大气廓线, 受其局限, 就不会是计算垂直气压梯度力精度高的非静力数值模式。

(5) 密度流试验表明, 在求得气块三维位移之后, 该三维位移隐含一个时间步长离散的、平均的三维散度场, 应以此绝热变率, 作用于气压、气温场, 本文密度流试验同样采用样条格式求取三维散度场。

(6) 密度流试验表明, 模式大气顶、底层气压场可取为静力平衡边界(符合大气运动物理原理), 其它变量场(暂)取为前/后差分边界, 其它数学与物理相结合的样条边界(如自然三次样条边界), 有待将来研究。

(7) 非静力全可压(完全弹性)密度流试验还表明, 密度流为"声波+重力波", 则声波与重力波能否分离, 怎样分离, 有待将来研究。

(8) 密度流试验还表明, 样条格式稳定性, 有赖于对变量场作适当的平滑, 但平滑是误差的重要来源, 本文实现了做"位温保守"气压、气温场平滑, 但如何依照三次样条拟合曲率判断做变量场(尤其风场)单点或局域平滑, 并且气压场平滑要保持质量守恒, 气温场平滑要保持能量守恒, 风场平滑要保持动量守恒, 但粘性物理作用除外, 因密度流试验是所谓"绝热无摩擦"理想试验, 则是平滑作用部分代替了粘性作用, 但在平滑之后, 容易求得平滑点/域的二阶可导补丁, 而应用三次样条光顺泛函数(代替平滑), 用于消除三次样条拟合"扰动"可能出现的不连续尖点或环绕有待将来研究。

(9) 密度流试验还表明, 样条格式动力框架误差来源于样条格式二阶空间余差和上游点路径达不到精确轨迹计算误差, 另外, 还存在确定模式大气顶层压、温场物理误差, 例如, 模式顶层气压场也是时变的, 其与大气受到非绝热加热有关(有日变化), 则如何确定模式顶层压、温场及其变化, 有待将来研究。