高階精度有限差分方案下的非跳點網格試驗：基于淺水波方程

徐道生 陳德輝 吳凱昕

中國氣象局廣州熱帶海洋氣象研究所/區域數值預報重點實驗室，廣州 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模式的預報性能。