中国气象学会主办。
文章信息
- 冯业荣, 薛纪善, 陈德辉, 吴凯昕. 2020.
- FENG Yerong, XUE Jishan, CHEN Dehui, WU Kaixin. 2020.
- GRAPES区域扰动预报模式动力框架设计及检验
- The dynamical core for GRAPES regional perturbation forecast model and verification
- 气象学报, 78(5): 805-815.
- Acta Meteorologica Sinica, 78(5): 805-815.
- http://dx.doi.org/10.11676/qxxb2020.048
文章历史
-
2019-05-15 收稿
2020-04-09 改回
2. 中国气象科学研究院,北京,100081;
3. 国家气象中心,北京,100081
2. Chinese Academy of Meteorological Sciences, Beijing 100081, China;
3. National Meteorological Center, Beijing 100081, China
增量形式的四维资料变分同化(4D-Var)系统是依据同化时间窗内的背景场(由非线性模式提供)和观测信息求代价函数的极小值,得到最优分析场。这个最优化过程需要反复积分切线性模式(TLM,Tangent Linear Model)及对应的伴随模式(ADM,Adjoint Model)。切线性模式是对非线性模式(NLM,Non-Linear Model)的线性化,国际上普遍采用代码对代码(Code to Code)的方式建立,即直接针对非线性模式代码进行线性化编程。这样产生的切线性模式与原来的非线性模式高度关联,在扰动很小的情形下具有很高的逼近精度。
另一种做法是构造扰动预报模式(PFM,Perturbation Forecast Model),该模式从组成非线性模式的原始微分方程组出发,通过小扰动线性化导出连续微分形式的扰动方程组,再采用一定的数学方案进行差分计算。研究(Lawless, et al, 2003)表明,尽管从扰动逼近无限小的角度来看,切线性模式比扰动模式具有更高的精度,但由于扰动量往往是有限大小的,故扰动模式能够代替切线性模式而得到较好的数学结果。Lawless等(2005)进一步试验表明,在增量4D-Var系统中使用扰动模式的收敛速率和分析结果与使用切线性模式几乎相同。
传统切线性模式对原来的非线性模式依赖性大,不具备相对独立性。非线性模式的更新或升级会导致切线性模式以及相应的伴随模式的重新编写,工作量巨大,维护成本较高。另外,切线性模式需要的内存空间也比较大,对于非线性模式中存在迭代计算的情形处理起来比较麻烦。而扰动模式独立性较好,不需要严格遵循复杂的非线性模式代码,直接基于非线性模式的基本物理原理和方案进行编写,能在代码层面保持清晰的物理意义,开发和维护成本较低。扰动模式另一个优点是可以避开非线性模式复杂计算方案带来的问题,例如针对非线性模式中半拉格朗日平流插值或某些迭代求解方案的切线性模式代码编写难度大,而扰动模式处理起来要容易得多。
扰动模式在研究和业务中已经得到较好的应用。英国气象局(UK Met Office)采用扰动模式建立了4D-Var系统(Lorenc, 1997;Lorenc, et al, 2005),并将该4D-Var系统与1.5 km水平分辨率的UKV模式结合,发展了一个主要用于预报对流降水的0—24 h短时临近预报系统,并应用于2012年伦敦奥运会,使得雷达、卫星、风廓线仪等高频次探测资料得到充分利用,降水预报效果明显优于基于3D-Var系统的区域模式和基于外推的临近预报系统(Ballard, et al,2016)。扰动模式还有其他的用途,例如通过扰动模式及其伴随模式评估模式预报对于观测资料的敏感性(Lorenc, et al,2014)。此外,还可用于分析模式误差的远距离传播特征等。
现阶段GRAPES区域模式的区域同化系统仍主要采用3D-Var,一个主要原因是还没有建立起GRAPES区域模式的较成熟的切线性模式以及相应的伴随模式。针对这一问题,不同于GRAPES全球模式的切线性模式编写方法,文中针对GRAPES区域模式采用小扰动方法开发线性模式动力框架,不仅为完整线性模式及其伴随模式的建立打下重要基础,而且由于该线性模式框架有很高的相对独立性,也可为进一步建立区域4D-Var系统提供便利。
文中从GRAPES模式的原始方程组(薛纪善等,2008)出发,推导线性扰动基本方程组,进而建立区域扰动预报模式— —GRAPES_PF。为评估GRAPES_PF对于GRAPES非线性模式的线性近似程度,设计了数值试验。通过同时在非线性模式和线性模式中施加一个中尺度初始扰动高压,分别模拟扰动的传播和演变。数值试验表明,GRAPES_PF较好地模拟了初始高压扰动所激发的惯性重力内波快速传播过程,且具有较高的计算精度,可作为GRAPES区域模式较好的线性模式。
2 扰动模式的设计 2.1 GRAPES非线性模式基本方程组GRAPES非线性模式在地形追随垂直坐标下的球面基本方程组如下(薛纪善等,2008):
运动方程
$ \begin{split} \frac{{\rm d}{{{V}}}}{{\rm d}t}=&-{c}_{p}\theta {\left({\nabla }_{{\rm{h}}}\varPi \right)}_{\hat {\textit z}}-{c}_{p}\theta \frac{\partial \varPi }{\partial \hat{{\textit z}}}\left({Z}_{{{s}}{{x}}}{{\bf{i}}}+{Z}_{{{s}}{{y}}}{{\bf{j}}}+{Z}_{{{s}}\hat{{{\textit z}}}}{{\bf{k}}}\right)+\\&\left\{fv+{\delta }_{{\rm{M}}}\left(\frac{vu\tan\varphi }{r}-\frac{uw}{r}\right)-{\delta }_{{\rm{\varphi}}}{f}_{{\rm{\varphi}}}w\right\}{{\bf{i}}}- \\ &\left\{fu+{\delta }_{{\rm{M}}}\left(\frac{{u}^{2}\tan\varphi }{r}+\frac{vw}{r}\right)\right\}{{\bf{j}}}-\\&\left\{g-{\delta }_{\rm M}\frac{{u}^{2}+{v}^{2}}{r}-{\delta }_{{\rm{\varphi}}}{f}_{{\rm{\varphi}}}u\right\}{{\bf{k}}}+{{{F}}}\\[-12pt] \end{split} $ | (1) |
式中,
连续方程
$ \frac{{{\rm{d}}\varPi }}{{{\rm{d}}t}} = - \hat \gamma \varPi {D_3} + \frac{{\hat \gamma }}{{{c_{{p}}}\theta }}\left({{Q_{\rm{R}}} + L{P_{\rm{w}}}} \right) $ | (2) |
式中,
散度
热力学方程
$ \frac{{\rm{d}}\theta }{{\rm{d}}t}=\frac{{Q}_{{\rm{R}}}+L{P}_{{\rm{w}}}}{{c}_{{{p}}}\varPi } $ | (3) |
水汽方程
$ \frac{{\rm{d}}q}{{\rm{d}}t}=-{P}_{{\rm{w}}} $ | (4) |
设预报变量
$\begin{split} {\left(\frac{{\rm{d}}A}{{\rm{d}}t}\right)}'&={\left(\frac{\partial A}{\partial t}+{{{V}}}{\text •} \nabla A\right)}'\approx \frac{\partial {A}'}{\partial t}+{{{V}}}{\text •} \nabla {A}'+{{{{V}}}}'{\text •} \nabla A\\&=\frac{{\rm{D}}{A}'}{{\rm{D}}t}+{{{{V}}}}'{\text •} \nabla A \end{split}$ | (5) |
上式表明,时间全导数的扰动等于扰动量的时间全导数加上扰动风场对基本态的平流。为简洁,推导过程已忽略基本态的标记“—”。其中全导数算子:
将式(1)写成分量形式并线性化
$ \frac{{{\rm{D}}u}'}{{\rm{D}}t}=-{{{{V}}}'}{\text •} \nabla u+{c}_{{{p}}}{{L}_{{\text{π}}{{x}}}\theta }'+{L}_{{\rm{\theta}}}({\varPi }_{x}'+{Z}_{sx}{\varPi }_{\hat{{\textit z}}}')+fv' $ | (6) |
$ \frac{{{\rm{D}}v}'}{{\rm{D}}t}={-{{{{V}}}'}{\text •} \nabla v+c}_{{{p}}}{L}_{{\text{π}}{y}}{\theta }'+{L}_{{\rm{\theta}}}({\varPi }_{y}'+{Z}_{{sy}}{\varPi }_{\hat{{\textit z}}}')-fu' $ | (7) |
$ \frac{{\rm{D}}{w}'}{{\rm{D}}t}=-{{{{V}}}'}{\text •} \nabla w+{c}_{{{p}}}{L}_{{\text{π}}\hat{{{\textit z}}}}\theta '+{L}_{{\rm{\theta}}}{\textit Z}_{{{s}}\hat{{{\textit z}}}}{\varPi }_{\hat{{\textit z}}}' $ | (8) |
式中,算子
对式(2)进行线性化,得到扰动连续方程
$ \frac{{\rm{D}}{\varPi }'}{{\rm{D}}t}=-{{{{V}}}'}{\text •} \nabla \varPi -{\alpha }_{{\rm{D}}}{\varPi }'-{\alpha }_{\Pi}{D}_{3}'+{\gamma }_{{\rm{\theta}}}{P'_{{\rm{w}}}} $ | (9) |
此处做简化,
对式(3)和(4)做线性化,分别得到扰动热力学方程和扰动水汽方程
$ \frac{{\rm{D}}{\theta }'}{{\rm{D}}t}=-{{{{V}}}'}{\text •} \nabla \theta +{\gamma }_{\Pi}{P'_{{{{\rm{w}}}}}} $ | (10) |
此处做简化
$ \frac{{{\rm{D}}q}'}{{\rm{D}}t}=-{{{{V}}}'}{\text •} \nabla q-{P'_{{\rm{w}}}} $ | (11) |
最终得到包含6个扰动变量(u′、v′、
$ \begin{split} \frac{{{\rm{D}}u}'}{{\rm{D}}t}=&-\left({u}'{u}_{x}+{v}'{u}_{y}+{\hat{w}}'{u}_{\hat{{\textit z}}}\right)+{c_p}{{L}_{{\text{π}}x}\theta }'+\\&{L}_{{\rm{\theta}}}({\varPi }_{x}'+{Z}_{sx}{\varPi }_{\hat{{\textit z}}}')+fv' \end{split}$ | (12) |
$\begin{split} \frac{{{\rm{D}}v}'}{{\rm{D}}t}=&{-\left({u}'{v}_{x}+{v}'{v}_{y}+{\hat{w}}'{v}_{\hat{{\textit z}}}\right)+c_p}{L}_{{\text{π}}{\rm{y}}}{\theta }'+\\&{L}_{{\rm{\theta}}}({\varPi }_{y}'+{Z}_{sy}{\varPi }_{\hat{{\textit z}}}')-fu' \end{split}$ | (13) |
$ \frac{{\rm{D}}{w}'}{{\rm{D}}t}=-\left({u}'{w}_{x}+{v}'{w}_{y}+{\hat{w}}'{w}_{\hat{{\textit z}}}\right)+{c_pL}_{{\text{π}}\hat{{{\textit z}}}}\theta '+{L}_{{{\theta}}}{Z}_{{{s}}\hat{{{\textit z}}}}{\varPi }_{\hat{{\textit z}}}' $ | (14) |
$\begin{split} \frac{{\rm{D}}{\varPi }'}{{\rm{D}}t}=&-\left({u}'{\varPi }_{x}+{v}'{\varPi }_{y}+{\hat{w}}'{\varPi }_{\hat{{\textit z}}}\right)-{\alpha }_{{\rm{D}}}{\varPi }'-\\&{\alpha }_{\Pi}{D}_{3}'+{\gamma }_{{\rm{\theta}}}{P}_{{\rm{w}}}'\end{split} $ | (15) |
$ \frac{{\rm{D}}{\theta }'}{{\rm{D}}t}=-\left({u}'{\theta }_{x}+{v}'{\theta }_{y}+{\hat{w}}'{\theta }_{\hat{{\textit z}}}\right)+{\gamma }_{\Pi}{P}_{{\rm{w}}}' $ | (16) |
$ \frac{{{\rm{D}}q}'}{{\rm{D}}t}=-\left({u}'{q}_{x}+{v}'{q}_{y}+{\hat{w}}'{q}_{\hat{{\textit z}}}\right)-{P}_{{\rm{w}}}' $ | (17) |
其中,
$ \hat w' = {Z_{sx}}u' + {Z_{sy}}v' + {Z_{{{s\hat {\textit z}}}}}w'\qquad\qquad\qquad $ |
与GRAPES非线性模式类似,采用两个时间层的半隐式半拉格朗日(Temperton, et al,1987))时间积分方案。
假定扰动量
$ \frac{{\left({A}'\right)}_{{\rm{a}}}^{{{n+1}}}-{\left({A}'\right)}_{{\rm{d}}}^{n}}{\Delta t}=\alpha {\left({L}_{{\rm{A}}}\right)}_{{\rm{a}}}^{n+1}+(1-\alpha) {\left({L}_{{\rm{A}}}\right)}_{{\rm{d}}}^{n} $ |
或写成
$ {\left({A}'\right)}_{{\rm{a}}}^{n+1}-{\left({A}'\right)}_{{\rm{d}}}^{n}={\alpha }_{{\rm{\varepsilon}}}{\left({L}_{{\rm{A}}}\right)}_{{\rm{a}}}^{{n}+1}+ {\beta }_{{\rm{\varepsilon}}}{\left({L}_{{\rm{A}}}\right)}_{{\rm{d}}}^{{n}} $ | (18) |
这里
因此,对扰动纬向速度方程(12),利用式(18)得到
$ {\left({u}'\right)}_{{\rm{a}}}^{{{n}}+1}-{\alpha }_{{\rm{\varepsilon}}}{\left[{L}_{{\rm{\theta}}}{(\varPi }_{x}'+{Z}_{{{sx}}}{\varPi }_{\hat{{\textit z}}}')+fv'\right]}_{{\rm{a}}}^{{{n}}+1}={X}_{u}' $ | (19) |
为书写简便,下面除必须外,下标a和时间上标n+1一概略去。于是上式变为
$ {u}'-{\alpha }_{{\rm{\varepsilon}}}f{v}'-{\alpha }_{{\rm{\varepsilon}}}{L}_{{\rm{\theta}}}{(\varPi }_{x}'+{Z}_{{{sx}}}{\varPi }_{\hat{{\textit z}}}')={X}_{{{u}}}' $ | (20) |
右端项
$\begin{split} {X}_{u}'=&({u}')_{{\rm{d}}}^{{{n}}}+{\beta }_{{\rm{\varepsilon}}}f{\left({v}'\right)}_{{\rm{d}}}^{{{n}}}+{\beta }_{{\rm{\varepsilon}}}{L}_{{\rm{\theta}}}{{(\varPi }_{x}'+{Z}_{{{s}}{{x}}}{\varPi }_{\hat{{\textit z}}}')}_{{\rm{d}}}^{{{n}}}+\\&\Delta t{c}_{{{p}}}{L}_{{\text{π}}{{x}}}{\theta }'^{{{n}}}-\Delta t{\left({{{V}}'}{\text •} \nabla u\right)}^{{{n}}}\end{split} $ |
类似地,扰动经向速度方程为
$ {v}'+{\alpha }_{{\rm{\varepsilon}}}f{u}'-{\alpha }_{{\rm{\varepsilon}}}{L}_{{\rm{\theta}}}{(\varPi }_{y}'+{Z}_{{\rm{s}}{\rm{y}}}{\varPi }_{\hat{{\textit z}}}')={X}_{v}' $ | (21) |
式中,
扰动垂直速度方程
$ {w}'+{\alpha }_{{\rm{\varepsilon}}}{c}_{{{p}}}\theta {'Z}_{{{s}}\hat{{{\textit z}}}}{\varPi }_{\hat{{\textit z}}}+{\alpha }_{{\rm{\varepsilon}}}{c}_{{{p}}}\theta {Z}_{{s}\hat{{\textit z}}}{\varPi }_{\hat{{\textit z}}}'={X}_{{{w}}}' $ | (22) |
式中,
将地形追随坐标垂直速度
$ \begin{split}&\hat w' - {Z_{{{sx}}}}u' - {Z_{{{sy}}}}v' + {\alpha _\varepsilon }{c_{{p}}}\theta {Z_{{{s\hat z}}}^2}\varPi {'_{\hat {\textit z}}} + {\alpha _\varepsilon }{c_{{p}}}{Z_{{{s\hat z}}}^2}{\varPi _{\hat {\textit z}}}\theta ' = X_{{{\hat w}}}' \end{split}$ | (23) |
式中,
扰动连续方程
$ \begin{split}& {\varPi }'+{\alpha }_{{\rm{\varepsilon}}}\left({u}'{\varPi }_{x}+{v}'{\varPi }_{y}+{\hat{w}}'{\varPi }_{\hat{\textit z}}\right)+{\alpha }_{{\rm{\varepsilon}}}{\alpha }_{{\rm{D}}}{\varPi }'+{\alpha }_{{\rm{\varepsilon}}}{\alpha }_{\Pi}{D}_{3}'={X}_{\varPi }' \end{split}$ | (24) |
式中,
潜热加热项
将
$\begin{split}& {\varPi }'={A}_{\varPi1}u'+{A}_{\varPi2}v'+{A}_{\varPi3}{\hat{w}}'+\\&{A}_{\varPi4}\left({u}'_{x}+{v}'_{y}+\hat{w}'_{\hat{{\textit z}}}\right)+{A}_{\varPi0} \end{split}$ | (25) |
其中系数
$\begin{split} &{A_{{{\varPi }}1}} = {\alpha _{\rm{\varepsilon }}}\dfrac{{{\alpha _{{\varPi }}}\dfrac{{\partial {Z_{sx}}}}{{\partial \hat {\textit z}}} - {\varPi _x}}}{{1 + {\alpha _{\rm{\varepsilon }}}{\alpha _{\rm{D}}}}},\quad{A_{{{\varPi }}2}} = {\alpha _\varepsilon }\dfrac{{{\alpha _{{\varPi }}}\dfrac{{\partial {Z_{sy}}}}{{\partial \hat {\textit z}}} - {\varPi _y}}}{{1 + {\alpha _{\rm{\varepsilon }}}{\alpha _{\rm{D}}}}},\\&{A_{{{\varPi }}3}} = - \frac{{{\alpha _\varepsilon }{\varPi _{\hat {\textit z}}}}}{{1 + {\alpha _{\rm{\varepsilon }}}{\alpha _{\rm{D}}}}},\quad{A_{{{\varPi }}4}} = - \frac{{{\alpha _{\rm{\varepsilon }}}{\alpha _{{\varPi }}}}}{{1 + {\alpha _{\rm{\varepsilon }}}{\alpha _{\rm{D}}}}},\quad{A_{{{\varPi }}0}} = \frac{{X_{{\varPi }}'}}{{1 + {\alpha _{\rm{\varepsilon }}}{\alpha _{\rm{D}}}}} \end{split}$ |
扰动水汽方程
$ {q}'=-{\alpha }_{{\rm{\varepsilon}}}{q}_{\hat{{\textit z}}}{\hat{w}}'+{X}_{{\rm{q}}}' $ | (26) |
式中,
这里扰动风对基态水汽的垂直平流用半隐式计算,水平平流用显式计算。
扰动位温方程形式与水汽方程完全类似,即
$ {\theta }'=-{\alpha }_{{\rm{\varepsilon}}}{\hat{w}}'{\theta }_{\hat{{\textit z}}}+{X}_{{\rm{\theta}}}' $ | (27) |
式中,
通过一系列代数消元后,扰动方程组最终化为只含变量
$ u'={A}_{{u}1}{\varPi }_{x}'+{A}_{{u}2}{\varPi }_{y}'+{A}_{{u}3}{\varPi }_{\hat{{\textit z}}}'+{A}_{{u}0} $ | (28) |
$ v'={A}_{{v}1}{\varPi }_{x}'+{A}_{{v}2}{\varPi }_{y}'+{A}_{{v}3}{\varPi }_{\hat{{\textit z}}}'+{A}_{{v}0} $ | (29) |
$\begin{split} \text{其中,} &{A_{{{u}}1}} = \frac{{{\alpha _{\rm{\varepsilon }}}{L_{\rm{\theta }}}}}{{1 + {{({\alpha _{\rm{\varepsilon }}}f)}^2}}},\quad{A_{{{u}}2}} = \frac{{{\alpha _{\rm{\varepsilon }}}^2f{L_{\rm{\theta }}}}}{{1 + {{({\alpha _{\rm{\varepsilon }}}f)}^2}}},\\&{A_{{{u}}3}} = \frac{{{\alpha _{\rm{\varepsilon }}}{L_{\rm{\theta }}}\left({{Z_{{sx}}} + {\alpha _{\rm{\varepsilon }}}f{Z_{sy}}} \right)}}{{1 + {{({\alpha _{\rm{\varepsilon }}}f)}^2}}},\quad{A_{{{u}}0}} = \frac{{{\alpha _\varepsilon }fX_{{v}}' + X_{{u}}'}}{{1 + {{({\alpha _{\rm{\varepsilon }}}f)}^2}}}\end{split} $ |
$\begin{split} &{A_{{{v}}1}} = - \frac{{{\alpha _{\rm{\varepsilon }}}^2f{L_{\rm{\theta }}}}}{{1 + {{({\alpha _{\rm{\varepsilon }}}f)}^2}}},\quad{A_{{{v}}2}} = \frac{{{\alpha _{\rm{\varepsilon }}}{L_{\rm{\theta }}}}}{{1 + {{({\alpha _{\rm{\varepsilon }}}f)}^2}}},\\ &{A_{{{v}}3}} = \frac{{{\alpha _{\rm{\varepsilon }}}{L_{\rm{\theta }}}\left({{Z_{{{sy}}}} - {Z_{{sx}}}{\alpha _{\rm{\varepsilon }}}f} \right)}}{{1 + {{({\alpha _{\rm{\varepsilon }}}f)}^2}}},\quad{A_{{{v}}0}} = \frac{{X_v' - {\alpha _\epsilon }fX_u'}}{{1 + {{({\alpha _{\rm{\varepsilon }}}f)}^2}}}\end{split} $ |
联立式(23)和(27),得
$ \begin{split} {\hat{w}}'=&\frac{{Z}_{{{sx}}}}{1-{\alpha }_{\varepsilon }^{2}{c}_{{{p}}}{Z}_{{{s}}\hat{{{\textit z}}}}^{2}{\varPi }_{\hat{{\textit z}}}{\theta }_{\hat{{\textit z}}}}{u}'+\frac{{Z}_{{{sy}}}}{1-{\alpha }_{\varepsilon }^{2}{c}_{{{p}}}{Z}_{{s}\hat{{{\textit z}}}}^{2}{\varPi }_{\hat{{\textit z}}}{\theta }_{\hat{{\textit z}}}}{v}'-\\ &\frac{{\alpha }_{\varepsilon }{c}_{{{p}}}\theta {Z}_{{s}\hat{{\textit z}}}^{2}}{1-{\alpha }_{\varepsilon }^{2}{c}_{{{p}}}{Z}_{{s}\hat{{\textit z}}}^{2}{\varPi }_{\hat{{\textit z}}}{\theta }_{\hat{{\textit z}}}}{{\varPi }}'_{\hat{{\textit z}}}+\frac{{X}_{\hat{{{w}}}}'-{\alpha }_{{\rm{\varepsilon}}}{c}_{{{p}}}{Z}_{{{s}}\hat{{{\textit z}}}}^{2}{\varPi }_{\hat{{\textit z}}}{X}_{\theta }'}{1-{\alpha }_{\varepsilon }^{2}{c}_{{{p}}}{Z}_{{{s}}\hat{{{\textit z}}}}^{2}{\varPi }_{\hat{{\textit z}}}{\theta }_{\hat{{\textit z}}}} \end{split}$ |
将式(28)和(29)代入上式,得到
$ {\hat{w}}'={A}_{\hat{{{w}}}1}{\varPi }_{x}'+{A}_{\hat{{{w}}}2}{\varPi }_{y}'+{A}_{\hat{{{w}}}3}{\varPi }_{\hat{{\textit z}}}'+{A}_{\hat{{{w}}}0} $ | (30) |
其中,
$ \begin{split} &\ {A_{{{\hat w}}2}} = \frac{{{Z_{{{sx}}}}{A_{{{u}}2}} + {Z_{{{sy}}}}{A_{{{v}}2}}}}{{1 - \alpha _\varepsilon ^2{c_{{p}}}Z_{{{s\hat {\textit z}}}}^2{\varPi _{\hat {\textit z}}}{\theta _{\hat {\textit z}}}}},\\ &\ {A_{{{\hat w}}3}} = \frac{{{Z_{{{sx}}}}{A_{{{u}}3}} + {Z_{{sy}}}{A_{{{v}}3}} - {\alpha _\varepsilon }{c_{{p}}}\theta {Z_{{{s\hat {\textit z}}}}^2}}}{{1 - \alpha _\varepsilon ^2{c_{{p}}}{Z_{{s\hat {\textit z}}}^2}{\varPi _{\hat {\textit z}}}{\theta _{\hat {\textit z}}}}},\\ &\ {A_{{{\hat w}}0}} = \frac{{{Z_{{{sx}}}}{A_{{{u}}0}} + {Z_{{{sy}}}}{A_{{{v}}0}} + X_{{{\hat w}}}' - {\alpha _{\rm{\varepsilon }}}{c_{{p}}}{Z_{{s\hat {\textit z}}}^2}{\varPi _{\hat {\textit z}}}X_\theta '}}{{1 - \alpha _\varepsilon ^2{c_{{p}}}{Z_{{{s\hat {\textit z}}}}^2}{\varPi _{\hat {\textit z}}}{\theta _{\hat {\textit z}}}}} \end{split} $ |
进一步将式(30)代入式(26)和(27)得到
$ {\theta }'={A}_{{{\theta}}1}{\varPi }_{x}'+{A}_{{{\theta}}2}{\varPi }_{y}'+{A}_{{{\theta}}3}{\varPi }_{\hat{{\textit z}}}'+{A}_{{{\theta}}0} $ | (31) |
$ \begin{split} \text{其中,}\qquad &{A_{{{\theta }}1}} = - {\alpha _{{\varepsilon }}}{\theta _{{{\hat {\textit z}}}}}{A_{{{\hat w}}1}},\quad{A_{{{\theta }}2}} = - {\alpha _{{\varepsilon }}}{\theta _{{{\hat {\textit z}}}}}{A_{{{\hat w}}2}},\qquad\qquad\\ &{A_{{{\theta }}3}} = - {\alpha _{{\varepsilon }}}{\theta _{{{\hat {\textit z}}}}}{A_{{{\hat w}}3}},\quad{A_{{\rm{\theta }}0}} = - {\alpha _{\rm{\varepsilon }}}{\theta _{{{\hat {\textit z}}}}}{A_{{{\hat w}}0}} + X_{{\theta }}'\\ &q' = {A_{{{q}}1}}\varPi _x' + {A_{{{q}}2}}\varPi _y' + {A_{{{q}}3}}\varPi _{\hat {\textit z}}' + {A_{{{q}}0}}\\[-10pt] \end{split} $ | (32) |
$ \begin{split} \text{其中,}\qquad &{A_{{{q}}1}} = - {\alpha _{{\varepsilon }}}{q_{{{\hat {\textit z}}}}}{A_{{{\hat w}}1}},\quad{A_{{{q}}2}} = - {\alpha _{{\varepsilon }}}{q_{{{\hat {\textit z}}}}}{A_{{{\hat w}}2}},\\ &{A_{{{q}}3}} = - {\alpha _{{\varepsilon }}}{q_{{{\hat {\textit z}}}}}{A_{{{\hat w}}3}},\quad{A_{q0}} = - {\alpha _\varepsilon }{q_{{{\hat {\textit z}}}}}{A_{{{\hat w}}0}} + X_{{q}}'\end{split} \hspace{27pt} $ |
方程(28)—(32)中,扰动变量(u'、v'、
沿用GRAPES非线性模式的空间离散差分方案,垂直方向采用Charney-Phillips跳点,水平方向采用Arakawa C网格。为提高精度,垂直差分使用非均匀分层的二阶精度差分方案。
空间离散化后的运动方程为
$ {u}'={A}_{{{u}}1}{\varPi }_{x}'+{A}_{{{u}}2}{\left({\overline{{\varPi }'}}^{yx}\right)}_{y}+{A}_{{{u}}3}{\left({\overline{{\varPi }'}}^{\hat{{\textit z}}x}\right)}_{\hat{{\textit z}}}+{A}_{{{u}}0} $ |
$ {v}'={A}_{{{v}}1}{\left({\overline{{\varPi }'}}^{xy}\right)}_{x}+{A}_{{{v}}2}{\varPi }_{y}'+{A}_{{{v}}3}{\left({\overline{{\varPi }'}}^{y\hat{{\textit z}}}\right)}_{\hat{{\textit z}}}+{A}_{{{v}}0} $ |
$ {\hat{w}}'={A}_{\hat{{{w}}}1}{\left({\overline{{\varPi }'}}^{x\hat{{\textit z}}}\right)}_{x}+{A}_{\hat{{{w}}}2}{\left({\overline{{\varPi }'}}^{y\hat{{\textit z}}}\right)}_{y}+{A}_{\hat{{{w}}}3}{\varPi }_{\hat{{\textit z}}}'+{A}_{\hat{{{w}}}0} $ |
变量上方的横线表示因为跳点的原因需要在相应的平面上取平均。将上面3个公式代入式(25)并离散化,得
$ \begin{split} {\varPi }'=&{A}_{{\varPi}1}{\overline{\left[{A}_{{{u}}1}{\varPi }_{x}'+{A}_{{{u}}2}{\left({\overline{{\varPi }'}}^{yx}\right)}_{y}+{A}_{{{u}}3}{\left({\overline{{\varPi }'}}^{\hat{{\textit z}}x}\right)}_{\hat{{\textit z}}}+{A}_{u0}\right]}}^{x}+\\ &{A}_{\varPi2}{\overline{\left[{A}_{{{v}}1}{\left({\overline{{\varPi }'}}^{xy}\right)}_{x}+{A}_{{{v}}2}{\varPi }_{y}'+{A}_{{{v}}3}{\left({\overline{{\varPi }'}}^{y\hat{{\textit z}}}\right)}_{\hat{{\textit z}}}+{A}_{{{v}}0}\right]}}^{y}+\\ &{A}_{\varPi3}{\overline{\left[{A}_{\hat{{{w}}}1}{\left({\overline{{\varPi }'}}^{x\hat{{\textit z}}}\right)}_{x}+{A}_{\hat{{{w}}}2}{\left({\overline{{\varPi }'}}^{y\hat{{\textit z}}}\right)}_{y}+{A}_{\hat{{{w}}}3}{\varPi }_{\hat{{\textit z}}}'+{A}_{\hat{w}0}\right]}}^{\hat{{\textit z}}}+ \\&{A}_{\varPi4}\left\{{\left[{A}_{{{u}}1}{\varPi }_{x}'+{A}_{{{u}}2}{\left({\overline{{\varPi }'}}^{yx}\right)}_{y}+{A}_{{{u}}3}{\left({\overline{{\varPi }'}}^{\hat{{\textit z}}x}\right)}_{\hat{{\textit z}}}+{A}_{u0}\right]}_{x}\right.+ \\&{\left[{A}_{{{v}}1}{\left({\overline{{\varPi }'}}^{xy}\right)}_{x}+{A}_{{{v}}2}{\varPi }_{y}'+{A}_{{{v}}3}{\left({\overline{{\varPi }'}}^{y\hat{{\textit z}}}\right)}_{\hat{{\textit z}}}+{A}_{{{v}}0}\right]}_{y}+ \\&\left.{\left[{A}_{\hat{w}1}{\left({\overline{{\varPi }'}}^{x\hat{{\textit z}}}\right)}_{x}+{A}_{\hat{{\rm{w}}}2}{\left({\overline{{\varPi }'}}^{y\hat{{\textit z}}}\right)}_{y}+ {A}_{\hat{w}3}{\varPi }_{\hat{{\textit z}}}'+{A}_{\hat{{{w}}}0}\right]}_{\hat{{\textit z}}}\right\}+{A}_{\varPi 0} \end{split}$ | (33) |
将式(33)在网格点上展开,得到如下形式的亥姆霍兹方程
$ \begin{split} & {B}_{1}{\varPi }_{{{i}},{{j}},{{k}}}+{B}_{2}{\varPi }_{{{i}}-1,{{j}},{{k}}}+{B}_{3}{\varPi }_{{{i}}+1,{{j}},{{k}}}+{B}_{4}{\varPi }_{{{i}},{{j}}-1,{{k}}}+{B}_{5}{\varPi }_{{{i}},{{j}}+1,{{k}}}+\\ &\quad{B}_{6}{\varPi }_{{{i}}+1,{{j}}+1,{{k}}}+{B}_{7}{\varPi }_{{{i}}+1,{{j}}-1,{{k}}}+{B}_{8}{\varPi }_{{{i}}-1,{{j}}-1,{{k}}}+{B}_{9}{\varPi }_{{{i}}-1,{{j}}+1,{{k}}}+\\ &\quad{B}_{10}{\varPi }_{{{i}},{{j}},{{k}}-1}+{B}_{11}{\varPi }_{{{i}}-1,{{j}},{{k}}-1}+{B}_{12}{\varPi }_{{{i}}+1,{{j}},{{k}}-1}+{B}_{13}{\varPi }_{{{i}},{{j}}-1,{{k}}-1}+\\ &\quad{B}_{14}{\varPi }_{{{i}},{{j}}+1,{{k}}-1} +{B}_{15}{\varPi }_{{{i}},{{j}},{{k}}+1}+{B}_{16}{\varPi }_{{{i}}-1,{{j}},{{k}}+1}+{B}_{17}{\varPi }_{{{i}}+1,{{j}},{{k}}+1}+\\&\quad{B}_{18}{\varPi }_{{{i}},{{j}}-1,{{k}}+1}+{B}_{19}{\varPi }_{{{i}},{{j}}+1,{{k}}+1}={B}_{0}\\[-10pt] \end{split} $ | (34) |
式中系数
亥姆霍兹方程(34)的求解采用广义共轭余差法(GCR)迭代求解,关于该方法详见薛纪善等(2008)。当前求解过程取
为使方程(34)的求解闭合,使用了上下边界为刚体的边界条件,即
上边界
$ {\left({\hat{w}}'\right)}_{{\rm{T}}}={\left[{A}_{\hat{{{w}}}1}{\left({\overline{{\varPi }'}}^{x\hat{{\textit z}}}\right)}_{x}+{A}_{\hat{{{w}}}2}{\left({\overline{{\varPi }'}}^{y\hat{{\textit z}}}\right)}_{y}+{A}_{\hat{{{w}}}3}{\varPi }_{\hat{{\textit z}}}'+{A}_{\hat{{{w}}}0}\right]}_{{\rm{T}}}=0 $ | (35) |
下边界
$ {\left({\hat{w}}'\right)}_{{\rm{B}}}={\left[{A}_{\hat{{{w}}}1}{\left({\overline{{\varPi }'}}^{x\hat{{\textit z}}}\right)}_{x}+{A}_{\hat{{{w}}}2}{\left({\overline{{\varPi }'}}^{y\hat{{\textit z}}}\right)}_{y}+{A}_{\hat{{{w}}}3}{\varPi }_{\hat{{\textit z}}}'+{A}_{\hat{w}0}\right]}_{{\rm{B}}}=0 $ | (36) |
为了测试扰动模式动力框架的合理性,首先利用非线性模式产生检验所需要的扰动场。具体做法为,将非线性模式GRAPES的物理过程全部关闭(干模式),设计两组试验,分别以2015年5月16日00:00:00 UTC的ECMWF分析场作为初值进行积分。第一组试验(CTRL)直接以ECMWF分析场作为模式初始场进行时间积分;第二组试验(PERT)在ECMWF分析场上叠加一个中尺度小高压扰动后进行时间积分。试验PERT与试验CTRL的差值作为真值用于对扰动模式GRAPES_PF的积分结果进行检验。非线性模式GRAPES和扰动模式GRAPES_PF的水平分辨率均取0.12°(约12 km),垂直层为67层,模式顶高度为30 km,时间积分步长为60 s。模式范围为(16°—29.8ºN,96°—110.04ºE)(图1),侧边界采用Davies松弛边界处理方案(薛纪善等,2008),现阶段扰动模式在侧边界上简单取0,即令侧边界无扰动。
使用12 km分辨率的GRAPES-3DVar同化系统产生扰动高压。利用3D-Var中的单站探空分析功能,在模式面约6 km高度(约在模式25层)的(22.78°N,102.7°E)位置上进行气压要素的单点分析,得到水平方向尺度半径约120 km大致符合高斯分布的初始高压扰动(
扰动求解过程从n时刻初始扰动
可以预计,如果扰动模式设计合理,扰动的传播图像应与非线性模式两次试验的差值场相似。图2是模式积分11步后第20层的水平扰动风场,其中图2a是积分两组非线性模式的差值场,代表扰动的实际演变(视为真值);图2b是从扰动初值出发经过扰动模式积分11步后得到的结果,代表扰动模式的预报值。图2c是图2a与图2b的差值,代表扰动模式的预报与实际情况的差异。可以看到,扰动模式的预报与非线性模式的结果非常相似,差异较小,能够较好地模拟扰动风场的变化。
图3是非线性模式和扰动模式不同时刻在模式第25层的预报结果。可以看出,人为设置高压扰动后,大气运动的大尺度准平衡遭到破坏,初始没有风场支持的扰动高压迅速激发出了强的非地转风,以惯性重力内波的形式迅速向四周传播。这实际上是一个非常快的地转适应过程,叶笃正等(1965)以及曾庆存(1963)和曾庆存等(1980)从理论上对这一过程作了充分的研究。由图3可见,扰动模式非常逼真地模拟了这个气压场向风场适应的过程。伴随非地转扰动能量迅速向外弥散,扰动气压场迅速减弱,扰动风场也随之迅速趋于减弱,模式大气迅速恢复到原来的准平衡状态。
适应过程是准线性快过程,是以辐散辐合运动和垂直运动为特征的(曾庆存,1963),因此线性扰动模式必须能合理反映这个基本特征。从给出的积分第11步后6个扰动变量经过初始扰动高压中心所做的纬度高度剖面(图4)可见,初始高压扰动除了激发水平方向迅速传播的重力波涟漪,也引起了垂直方向的运动增量,并进而引起各变量在垂直方向上重新分配。与非线性模式相比,扰动模式能够较好地预报扰动随时间变化的空间分布,各个扰动变量的数量级和正负数值分布与非线性模式的结果相近。值得一提的是,
为了整体评估扰动模式相对于非线性模式的线性近似程度,利用Frobonius范数对扰动模式进行正确性检验。若将初始基本态(
$ {\rm{NLM}}\left({x}_{0}+\gamma {{\rm{\delta}}x}_{0}\right)={\rm{NLM}}\left({x}_{0}\right)+{R}_{0\to i}\gamma {{\rm{\delta}}x}_{0}+O\left({\gamma }^{2}\right) $ | (37) |
设
$ \begin{split} F\left(\gamma \right)&=\frac{\left\|{\rm{NLM}}\left({x}_{0}+\gamma {{\rm{\delta}}x}_{0}\right)-{\rm{NLM}}\left({x}_{0}\right)\right\|}{\left\|\rm{PFM}\left(\gamma {{\rm{\delta}}x}_{0}\right)\right\|} \\& = \frac{\left\|{R}_{0\to i}{{\rm{\delta}}x}_{0}\right\|}{\left\|{M}_{0\to i}{{\rm{\delta}}x}_{0}\right\|}+O\left(\gamma \right) \end{split} $ | (38) |
式中,‖•‖表示Frobonius范数。当
根据上、下边界条件公式(35)和(36),求解亥姆霍兹方程需要在模式顶和模式底计算扰动气压垂直导数
$ {\left({\varPi }_{\hat{{\textit z}}}'\right)}_{{\rm{bottom}}\;{\rm{or}}\;{\rm{top}}}=\frac{g{\theta }'}{{c}_{p}{Z}_{{{s}}\hat{{{ {\textit z}}}}}{\theta }^{2}} $ | (39) |
式中,
表1给出不同
参数γ | 虚层方案 | 静力平衡方案 |
1.0 | 0.8990750313 | 0.9092131257 |
0.1 | 1.038652658 | 0.9952505231 |
0.01 | 1.054983139 | 1.004807591 |
0.001 | 1.058226705 | 1.005833268 |
0.0001 | 1.192282557 | 1.037582755 |
为发展基于扰动预报模式的四维变分资料同化系统,从GRAPES数值预报模式的原始方程组出发,采用小扰动分离方法推导了连续微分形式的线性扰动预报方程组,开发了扰动预报模式GRAPES_PF动力框架。扰动模式中时间积分采用两个时间层的半隐式半拉格朗日方案,空间差分采用Arakawa C水平网格跳点和Charney/Phillips垂直跳层方案,利用代数消元法进一步推导,得到关于扰动无量纲气压(
为了测试所开发的扰动模式动力框架的合理性,在干模式(无物理过程)条件下开展了数值试验。首先在非线性模式中施加一个初始扰动,积分非线性模式得到扰动的时间演变场,然后将相同的初始扰动作为扰动模式的初值进行时间积分,最后比较扰动模式与非线性模式的预报结果。
试验在模式第25层(约6 km处)施加了一个水平半径约120 km的中尺度扰动高压,根据地转适应理论,在这样小的尺度下气压场将向风场适应,即小高压所激发的不稳定能量将迅速向四周频散,导致扰动趋于减弱。数值试验表明,初始的非地转平衡状态激发了一个快速向外传播的惯性重力内波,迅速引起了水平辐散风和铅直运动,进而引起位温和湿度扰动。随着时间的推移,适应过程使扰动气压场迅速频散消失,扰动风场也随之迅速减弱,大气运动重新恢复到原来的平衡状态。对比非线性模式和扰动模式,两者动力热力场分布特征相似,数值结果接近。地转适应过程是以水平辐散辐合运动和垂直运动为特征的准线性快过程,试验结果表明,扰动模式GRAPES_PF合理反映了这个运动的基本特性,可作为GRAPES非线性模式系统的近似线性模式,为下一步发展基于扰动模式的四维变分同化系统奠定了研究和开发基础。
文中仅设计开发了扰动模式的动力框架,完善该扰动模式还需要逐步完成各种线性物理过程方案的开发。
致 谢:袁丽君参与了扰动模式的检验以及部分代码编写和测试,吴亚丽进行了非线性模式的运算,在此一并感谢。
薛纪善, 陈德辉. 2008. 数值预报系统GRAPES的科学设计与应用. 北京: 科学出版社, 65-136. Xue J S, Chen D H. 2008. Scientific Design of the GRAPES Numerical Weather Prediction Model and Its Application. Beijing: Science Press, 65-136 (in Chinese)
|
叶笃正, 李麦村. 1965. 大气运动中的适应问题. 北京: 科学出版社, 107-124. Ye D Z, Li M C. 1965. On the Adaptation of the Atmospheric Motion. Beijing: Science Press, 107-124 (in Chinese)
|
曾庆存. 1963. 大气中的适应过程和发展过程(一): 物理分析和线性理论. 气象学报, 33(2): 163-174. Zeng Q C. 1963. The adaptation and evolution processes in atmosphere(I): Physical analysis and linear theory. Acta Meteor Sinica, 33(2): 163-174. (in Chinese) |
曾庆存. 1979. 数值天气预报的数学物理基础, 第一卷. 北京: 科学出版社, 457-464. Zeng Q C. 1979. The Physical-Mathematical Basis of Numerical Weather Prediction (I). Beijing: Science Press, 457-464 (in Chinese)
|
曾庆存, 叶笃正. 1980. 旋转大气中运动的适应过程. 力学学报, 16(1): 1-11. Zeng Q C, Ye D Z. 1980. The adaptation processes of motion in rotating atmosphere. Chinese J Theor Appl Mech, 16(1): 1-11. DOI:10.3321/j.issn:0459-1879.2005.01.001 (in Chinese) |
Ballard S P, Li Z H, Simonin D, et al. 2016. Performance of 4D-Var NWP-based nowcasting of precipitation at the Met Office for summer 2012. Quart J Roy Meteor Soc, 142(694): 472-487. DOI:10.1002/qj.2665 |
Lawless A S, Nichols N K, Ballard S P. 2003. A comparison of two methods for developing the linearization of a shallow-water model. Quart J Roy Meteor Soc, 129(589): 1237-1254. DOI:10.1256/qj.02.75 |
Lawless A S, Gratton S, Nichols N K. 2005. An investigation of incremental 4D-Var using non-tangent linear models. Quart J Roy Meteor Soc, 131(606): 459-476. DOI:10.1256/qj.04.20 |
Lorenc A C. 1997. Development of an operational variational assimilation scheme. J Meteor Soc Japan, 75(1): 339-346. |
Lorenc A C, Rawlins F. 2005. Why does 4D-Var beat 3D-Var?. Quart J Roy Meteor Soc, 131(613): 3247-3257. DOI:10.1256/qj.05.85 |
Lorenc A C, Marriott R T. 2014. Forecast sensitivity to observations in the Met Office Global numerical weather prediction system. Quart J Roy Meteor Soc, 140(678): 209-224. DOI:10.1002/qj.2122 |
Temperton C, Staniforth A. 1987. An efficient two-time-level semi-Lagrangian semi-implicit integration scheme. Quart J Roy Meteor Soc, 113(477): 1025-1039. DOI:10.1002/qj.49711347714 |