高阶精度有限差分方案下的非跳点网格试验：基于浅水波方程

徐道生 陈德辉 吴凯昕

中国气象局广州热带海洋气象研究所/区域数值预报重点实验室，广州 510640

1 引言

随着近年来各种观测资料的不断增加以及人们对大气变化规律的认识不断加深，数值天气预报模式的预报性能也得到了长足的进步。除了改进初始场和物理参数化方案之外，构造更加精确的动力框架（包括时间和空间离散化方案、网格设置、垂直分层设计等）也是提升模式预报效果的重要途径。

在原始方程大气模式中，地转适应过程的物理机制是重力惯性波对非地转能量的频散。因此模式对于地转适应过程的模拟能力，主要取决于它的差分格式能否正确地描写重力波的频散关系。大量研究结果表明（Winninghoff, 1968; Mesinger, 1973; 刘宇迪等, 2001; Rajpoot et al., 2012），网格点上变量分布对于地转适应过程的模拟近似程度有重要影响。Arakawa and Lanmb（1977）通过对几种不同变量分布的差分格式解进行分析，发现C网格在模拟地转适应过程方面的性能要明显优于A网格（即非跳点网格）、B网格、D网格和E网格。在变形半径能够被分辨的情况下，C网格可以比其他网格用更小的代价产生更优的数值离散，也就是说它可以用更优的性价比减少频散问题导致的计算噪音，目前C网格已经被广泛地应用于数值天气预报模式（Unified模式：Staniforth et al., 2006; GRAPES模式：薛纪善和陈德辉, 2008; WRF模式：Skamarock et al., 2008）。但是C网格也存在一些难以克服的缺陷：例如由于C网格需要将风场和其它模式变量错开分布，使得动力框架与物理过程的耦合变得更加复杂。因为在每次调用物理参数化方案之前需要将各个预报量插值到A网格，在完成物理过程反馈计算之后还要再次将它插值回到C网格。另外，C网格中水平风速分量u和v分布在不同网格上，在半隐式半拉格朗日框架下计算平流项和科氏力项的时候也需要进行多次插值计算。在C网格模式中这些反复插值的过程不仅增加了额外的计算量，也不可避免的在预报过程中引入了计算误差。

考虑到不同网格具有各自的独特优势，因此在很多模式中人们会基于不同的情况来选择其它几种网格（目前A、B、D、E网格仍然在模式中被采用），同时采取一些额外措施来缓解这些网格在频散计算方面的误差。Randall（1994）提出了一种利用涡度、散度和质量场作为预报变量的非跳点网格，通过和使用风场和质量场作为预报变量的C网格进行比较，发现新的非跳点网格对能量频散的模拟效果更优。宇如聪（1994）对E网格下的差分计算性质进行总结，并给出了一种解决E网格中两个C网格分离问题的办法。McGregor（2005）设计了一种非跳点网格，它仅仅需要在计算重力波的时候将风场可逆地插值到C网格上，通过浅水波方程试验发现这种网格的能量频散效果优于传统的C网格。Chen（2016）在涡度—散度和速度两种变量情况下对跳点网格与非跳点网格之间的等价关系进行了分析。Konor and Randall（2018）对各种水平和垂直网格的性能进行了全面的比较，认为Z网格和C网格更加适用于高分辨率的非静力模式。Miura（2019）将六边形C网格和改进的B网格进行比较，发现后者在纬向平衡流测试中具有更好的数值收敛速度，但是在其它方面则接近或差于C网格。Xie（2019）在Z网格的基础上提出一种广义Z网格，改进了Z网格变量离散化时的梯度等数值算法，使它在更大时间步长下也能保持计算稳定性。

近年来采用高阶精度算法对数值模式进行离散化已经成为主流趋势（Lin, 2004; Li et al., 2015）。已经有研究结果表明，在提高计算精度以后，非跳点网格的频散关系是有可能会发生变化的。Purser and Leslie（1988）基于高阶精度有限差分方案开展了非跳点网格的模拟试验，发现提高差分计算精度以后模式的模拟结果优于二阶精度跳点网格和低阶精度非跳点网格。最近Chen et al.（2018）则利用快速黎曼求解器构造了一个基于有限体积的动力框架，同样发现在高阶精度离散化方案下C网格与A网格的模拟结果差别并不明显（考虑粘性项的情况）。为了研究高阶精度方案下使用非跳点网格构造模式动力框架的可行性，本文基于浅水波方程在不同精度差分格式下使用跳点网格和非跳点网格进行了模拟试验。第2节首先给出跳点网格和非跳点网格下二阶、四阶和六阶精度的有限差分格式，第3节将这些不同精度的差分格式应用于浅水波方程的离散化，对相应的频散关系进行比较和分析。第4节针对虚假高频短波的处理方法进行介绍。在此基础上第5节利用一维浅水波方程开展数值模拟试验，最后在第6节进行总结和讨论。

2 高阶精度有限差分方案

传统差分格式得到某一点上高阶精度一阶导数表达式的方法是：首先确定达到要求精度时所需要的格点数，然后得出该点上含待定系数的导数表达式，最后利用Taylor展开法联立方程组求解待定系数。非跳点网格（Li, 2005）和跳点网格（Chu and Fan, 1997）下的二阶、四阶和六阶精度显式差分格式列于表1中。

3 两种网格的频散关系

地转适应过程基本上是一种线性过程，在不考虑地形影响和科氏力随纬度变化的情况下，正压大气原始方程组可以简化为如下形式的线性化二维浅水波方程：

按照Arakawa and Lamb（1977）的定义，在A网格中所有变量都分布在同一格点上，而C网格中u，v和z分别位于不同的格点上（如图2所示）。为了便于书写，定义如下舒曼算符：

分别将不同网格下各变量对应的半离散差分方程写为如下六种形式：

（1）A网格+二阶精度：

（2）A网格+四阶精度：

表1 跳点和非跳点网格下的不同精度差分格式Table 1 Difference schemes with different accuracies under staggered and nonstaggered grids

（3）A网格+六阶精度：

图1 二维浅水波方程中无量纲频率（等值线）和波数的真实关系，图中kd/pi和ld/pi分别表示x方向和y方向的无量纲波数（pi为圆周率）Fig. 1 Contours of the nondimensional frequency (contour) as a function of the nondimensional horizontal wavenumbers for a true solution of a shallow water equation. kd/pi and ld/pi are the nondimensional wavenumbers at xa nd y directions, respectively, pi is the ratio of the circumference of a circle to its diameter.

（4）C网格+二阶精度：

（5）C网格+四阶精度：

（6）C网格+六阶精度：

图2 浅水波试验中的二维网格变量分布：（a）A网格；（b）C网格Fig. 2 Distributions of variables in two-dimensional grid for the test of a shallow water equation: (a) A grid; (b) C grid

表2 A网格和C网格下不同精度差分格式的频率方程Table 2 Frequencies equation of Finite-difference grids A and C with different accuracies

总的来说，提高计算精度以后可以使得A网格在低波数区的模拟效果接近二阶精度C网格，但是不能避免频率极值的出现，只能使得极值的中心随着差分精度的提高向更高波数区域移动。这说明提高差分计算精度可以有效的减少低频波段（大于四倍格距）的频散误差，但是对于高频短波则没有改善。

4 高频短波的处理

从图3可以看出，在高阶差分精度下A网格的频率极值中心位于2～4倍格距波之间（即图3中 (l,k)d/(2pi)≤0.5的区域），也就是说它会导致高频短波的模拟出现能量无法频散的现象。由于大部分数值预报模式的实际有效分辨率都在6～8倍格距之间，为了避免“混淆误差”导致的非线性计算不稳定现象，通常需要对4倍格距以下的高频短波进行抑制。因此即使高阶精度差分算法下A网格仍然会出现高频极值，但是只要对它进行适当的处理，仍然可以对实际的波取得较好的模拟效果。例如Chen et al.（2018）在黎曼求解器中的隐式扩散项来处理高阶精度A网格的高频噪音之后，使得高阶精度非跳点网格的模拟效果达到和C网格相接近的程度。本文借鉴Chen et al.（2018）的做法，在预报方程中显式地增加一个粘性项来滤除高阶精度A网格引起的高频噪音。

图3 不同网格和差分精度下的无量纲频率（等值线）：（a）C网格＋二阶精度；（b）C网格＋四阶精度；（c）C网格＋六阶精度；（d）A网格＋二阶精度；（e）A网格＋四阶精度；（f）A网格＋六阶精度。深蓝色方框表示波长大于4倍格距的波数范围Fig. 3 Nondimensional frequency (Contours) under different grids and accuracy of differential computation : (a) C grid ＋ 2nd order; (b) C grid ＋ 4th order; (c) C grid ＋ 6th order; (d) A grid ＋ 2nd order; (e) A grid ＋ 4th order; (f) A grid ＋ 6th order. The blue square represents wavelength exceeding four-grid interval

5 数值试验

5.1 试验设计

为了对第3节的理论分析进行验证，下面通过如下形式的一维线性化浅水波方程（增加了粘性项）进行相应的数值模拟试验：

将初始扰动设置在-1000≤x≤1000 m范围内，并且是以上两个扰动的合成：

为了测试各种方案在不同波段下的模拟能力，分别选取十倍格距波和四倍格距波进行试验。各个对比试验的具体参数设置见表3。

表3 试验设计Table 3 Design of experiment

5.2 结果分析

在公式（12）至公式（15）中波数k分别取2π/(10Δx)和2π/(4Δx)时，即可得到波长等于十倍和四倍格距的初始扰动。图4为当t=300 s时不同精度下跳点网格和非跳点网格的十倍格距波模拟结果。对于完全可被分辨的十倍格距波来说，二阶精度C网格的模拟结果已经非常接近真实解（图4a），而二阶精度A网格的模拟结果会出现位相滞后现象（图4b），这是由于它的群速度偏慢导致的（Arakawa and Lamb, 1977）。但是当差分精度提高到四阶和六阶以后，A网格模拟结果（图4c、d）基本与二阶精度C网格模拟结果一致。这表明A网格下低频长波的模拟结果可以随计算精度的提高而得到明显的改善，并在四阶精度时达到与二阶精度C网格相接近的效果。这与图3中A网格低波数区的频率随计算精度提高而改进的结论是一致的。从均方根误差随时间的变化来看（图5），也可以看出高阶精度下A网格的模拟误差明显优于二阶精度，并与二阶精度C网格的模拟结果已经很接近。需要指出，模拟结果中出现的短波噪音主要由两种来源：（1）由于不连续性引起的噪音。它主要出现在两个波列之间的区域，这种噪音在A网格和C网格都会出现；（2）由于频散误差引起的噪音。它只出现在高阶精度A网格的模拟结果中，并且主要位于两列波的前方（见图4c和图4d中x=±5000 m处附近），它的频率会明显高于不连续性引起的噪音。通过在预报方程中显式增加粘性项以后，高阶精度差分方案下A网格模拟结果中的两种高频噪音都得到了较好地抑制（图5a）。从均方根误差来看（图5b），滤波以后A网格的误差（蓝色线）要小于滤波前（绿色线），并且更加接近二阶精度C网格的模拟精度（黑色线）。

图4 当t=300 s时在跳点网格和非跳点网格下的波长为十倍格距的重力波模拟：（a）test-C-2nd-10；（b）test-A-2nd-10；（c）test-A-4th-10；（d）test-A-6th-10Fig. 4 Simulation of gravity wave whose wavelength is set as 10 grid spacing using staggered and nonstaggered grids at t=300 s: (a) test-C-2nd-10;(b) test-A-2nd-10; (c) test-A-4th-10; (d) test-A-6th-10

图5 （a）当t=300 s时test-A-6th-10-f试验对重力波的模拟；（b）不同方案下波长为十倍格距的重力波均方根误差随时间的增长Fig. 5 (a) Simulation of gravity wave for test-A-6th-10-f at t=300 s; (b) the growth of the root-mean-squared error underdifferent schemes for the gravity wave whose wavelength is set as 10 grid spacing.

对于四倍格距波，二阶精度下C网格的模拟结果同样明显的优于二阶精度A网格（图6a，b）。由于在二阶精度下A网格的四倍格距波群速度正好为0，所以波一直在x=0附近做惯性震荡，基本无法模拟出它向两侧的能量频散过程。随着差分计算精度的提高，A网格对四倍格距波能量频散过程的模拟结果逐渐得到改进（图6c，d），这是由于计算精度提高以后A网格的频率极大值逐渐向更高波数区移动引起的。由于波的不连续性，所有模拟结果中两个列波之间的区域都出现了明显的高频噪音，而且高阶精度方案下A网格模拟结果中在波列的前方出现了更高频率的噪音（这主要是由于A网格的频散误差引起的），这与十倍格距波的模拟结果是一致的。对高频短波进行扩散以后（图7），可以看到几乎所有的信号都被耗散（包括四倍格距波本身和A网格频散误差导致的高频噪音），只剩下一些由于波动不连续性而导致的噪音，此时六阶精度A网格和二阶精度C网格的结果差别很小（图7a，b）。因此，虽然在高阶精度A网格对格点尺度短波的模拟存在较大误差，但是只要采用一定的特殊处理（比如通过增加粘性项进行滤除），就可以避免它对模式实际预报结果的影响。

从图3来看，在高频波段提高计算精度对减少C网格的频散误差也有一定的改进作用。为了测试这种改进对于重力波模拟结果的影响，分别采用四阶和六阶精度差分方案进行试验（图8）。和二阶精度模拟结果（图4a和图6a）进行比较，可以看出在C网格下高阶精度算法对模拟结果不但没有改进，反而会引入很多虚假的高频噪音（特别是在六阶精度的情况下）。在C网格下高阶精度算法虽然能够更准确地描述2～4倍格距波的频散关系，但是在积分过程中出现的这些格点尺度波通常属于不能被网格系统准确描述的高频噪音，它是造成非线性计算不稳定的主要原因。如果不对这些高频噪音进行耗散处理，高阶精度算法将会高估长波向短波的能量转换，导致短波能量的虚假增长。因此，在C网格下使用高阶精度差分算法对于实际波的模拟来说是没有意义的，这也是目前C网格模式中很少使用高阶精度有限差分算法的主要原因。

图6 当t=300 s时在跳点网格和非跳点网格下的波长为四倍格距的重力波模拟：（a）test-C-2nd-4；（b）test-A-2nd-4；（c）test-A-4th-4；（d）test-A-6th-4Fig. 6 Simulation of the gravity wave whose wavelength is set as 4 grid spacing using staggered grid and unstaggered at t=300s: (a) test-C-2nd-4;(b) test-A-2nd-4; (c) test-A-4th-4; (d) test-A-6th-4

图8 跳点网格下不同精度算法对重力波模拟：（a）test-C-4th-10；（b）test-C-6th-10；（c）test-C-4th-4；（d）test-C-6th-4Fig. 8 Simulation of gravity wave with different accuracies under staggered grids: (a) test-C-4th-10; (b) test-C-6th-10; (c) test-C-4th-4; (d) test-C-6th-4

6 总结

为了测试高阶精度有限差分格式下利用非跳点网格构建模式动力框架的可行性，本文基于浅水波方程，在高阶精度差分方案下对A网格和C网格进行分析和比较，并开展相应的模拟试验。主要得到以下结论：

（1）在二阶精度差分方案下，A网格的频散误差明显大于C网格。随着计算精度的提高，A网格出现频率极大值的位置逐渐向高频波段移动，低波数区的频散过程模拟结果更加接近真实解，非跳点网格与跳点网格的模拟结果差别也随之变小。在高波数区高阶精度A网格的频散误差仍然比较明显，而C网格的频散关系则随着计算精度的提高而进一步得到改进。

（2）一维浅水波数值试验结果表明：对于可以被模式完全识别的波（以十倍格距波为例）来说，在高阶精度算法下A网格的模拟结果可以接近二阶精度C网格的模拟结果。而对于格点尺度的高频短波（以四倍格距波为例），即使在六阶精度下A网格的模拟结果也仍然存在比较明显的误差，它对频散过程的模拟效果不如二阶精度C网格。

（3）通过显式增加粘性项对高频短波进行耗散以后，高阶精度A网格模拟结果中由于频散误差而出现的高频噪音得到了有效抑制，同时对四倍格距波的耗散程度也与二阶精度C网格的模拟结果相接近。

综上所述，通过使用粘性项对高频短波进行抑制的情况下，高阶精度A网格对重力波频散过程的模拟结果可以达到和二阶精度C网格相接近的模拟效果，因此在高阶精度差分格式下使用A网格构建模式动力框架是一种可行的做法。下一步我们将继续开展基于GRAPES模式的高阶精度差分和非跳点网格试验，进一步通过实际天气过程的预报来对本文的研究结论加以验证，从而最终改善GRAPES模式的预报性能。