1 引言

大气中包含多种波动，其中传播速度最快的是声波。对于预报方程采用显式时间差分方案的大气数值模式而言，模式的最大积分时间步长会受到传播速度最快的波动——声波的限制（Durran and Blossey, 2012）。由于声波没有显著的气象意义，为了避免声波对积分时步的限制，部分数值模式会采用滤除声波的控制方程组。在大气控制方程中滤除声波是指在一定的假设条件下，对大气控制方程做一定的简化和近似，滤除声波并保留具有气象意义的波动。模式中常用的滤除声波的控制方程有包辛尼斯克方程、滞弹性方程（Ogura and Phillips, 1962; Wilhelmson and Ogura, 1972; Lipps and Hemler, 1982; Bannon, 1996）和假不可压方程（Durran, 1989, 2008）。

包辛尼斯克近似忽略了连续方程中的密度个别变化项，从而使连续方程简化为了风场无辐散形式的诊断方程。包辛尼斯克近似仅适用于垂直尺度小于大气标高的大气运动，在描述深对流时，由于未考虑层结密度随高度的变化，则会产生较大误差。为了解决上述问题，滞弹性近似引入了一种新形式的连续方程，适用于研究深对流运动。

滞弹性方程组的构建可以追溯到1953年，Batchelor（1953）基于“密度和气压的分布与绝热分层的大气中的密度和气压分布接近”的假设给出了一组简化的控制方程，这就是最早的滞弹性方程组。Ogura and Phillips（1962）通过严格的尺度分析也得到了同样形式的大气控制方程组，并命名为“滞弹性”方程组。Dutton and Fichtl（1969）拓展了包辛尼斯克近似，使其适用于研究可压缩流体中的深对流运动。Gough（1969）在滞弹性方程中考虑了分子输送和辐射传输过程。Wilhelmson and Ogura（1972）在Ogura and Phillips（1962）工作的基础上详细分析了气压扰动对饱和水汽压的影响。Lipps and Hemler（1982）假设位温的基本态随高度缓慢变化，利用尺度分析得到了描述湿深对流的滞弹性方程组。随后，Lipps（1990）改进了Lipps and Hemler（1982）推导滞弹性方程组的过程，不再参考两个经验参数，使结果更加客观可信。Bannon（1996）建立了包含水汽和相变过程的能量守恒的湿滞弹性方程组。

Durran（1989）提出了一种新的滤除声波的动力学方程组，命名为“假不可压方程”，并成功用于构建数值模式。随后，Almgren et al.（2008）和O’Neill and Klein（2014）拓展了Durran（1989）的工作，前者建立了包含源汇项的假不可压模式，后者建立了包含水物质的假不可压模式。假不可压近似成立的条件有两个：一是扰动的拉格朗日时间尺度远大于声波传播的时间尺度，二是扰动气压远小于平均状态的气压。从物理意义上来看，假不可压近似的质量平衡中部分考虑了密度扰动的影响，即考虑与温度场扰动相关的密度扰动，而忽略与气压扰动相关的密度扰动。而在滞弹性近似中，密度扰动对质量平衡的影响则被完全忽略了。Smolarkiewicz and Dörnbrack（2008）以及Achatz et al.（2010）认为假不可压近似的动力场比滞弹性近似的动力场更接近完全可压缩大气。Nance and Durran（1994）认为假不可压系统的准确性高于滞弹性系统，Klein（2009）也证实了这一点。由此可见，假不可压近似比滞弹性近似更接近真实大气。此外，从成立条件来看，假不可压方程不涉及对扰动位温的量级的假设，而滞弹性方程需要，这使得假不可压近似的适用范围更广。综上所述，假不可压方程比滞弹性方程的准确性更高，适用范围更广，更适合用于构建数值预报模式。

WRF（Weather Research and Forecasting）模式广泛用于数值预报和研究，拥有成熟的平流方案和微物理模块。WRF模式ARW（Advanced Research WRF）动力框架采用了地形追随坐标系（η坐标系），可以简化复杂地形的下边界条件，能够较好地处理复杂地形问题。本文选取WRF-ARW动力框架作为模式的框架基础，用假不可压近似（Durran, 1989）改写模式的控制方程组，建立了基于ARW动力框架的假不可压模式。

本文将给出假不可压模式的建立过程，之后采用相同的试验设计方案，用假不可压模式与WRF模式，分别进行湿热泡对流试验和重力流流理想试验，并对两种模式的模拟结果进行分析和比较。

2 假不可压控制方程的建立

Durran提出的“假不可压方程”与完全可压的方程的核心区别在于连续方程的形式不同。在“假不可压方程”中，连续方程的形式如下：

$ \frac{{{\rm{d}}{\rho ^ * }}}{{{\rm{d}}t}} + {\rho ^ * }\left({\nabla \cdot {\mathit{\boldsymbol{V}}_z}} \right) = 0, $

(1)

其中，V_z是三维风矢量，ρ^*是假不可压密度，定义是${\rho ^ * } = (\bar \theta {\rm{/}}\theta)\bar \rho $或者$\bar \pi = {\left[ {({R_d}{\rm{/}}{p_0}){\rho ^ * }\theta } \right]^{\frac{{{R_d}}}{{{c_v}}}}}$。θ是位温。$\bar \pi $是平均无量纲气压，定义是$\bar \pi = {\left({\bar p{\rm{/}}{p_0}} \right)^{\frac{{{R_d}}}{{{c_p}}}}}$。p、ρ、θ是平均状态的气压、密度和位温，仅随高度变化。p₀是参考气压，取1000 hPa。R_d是干空气的比气体常数，c_p和c_v分别是干空气的等压和等容比热容。

2.1 连续方程的推导

由Durran提出的假不可压方程（1），利用z坐标系和η坐标系的转换关系（Kasahara, 1974）：

$ \frac{\partial }{{\partial z}} = \left({\frac{{\partial \eta }}{{\partial z}}} \right)\frac{\partial }{{\partial \eta }}, $

(2)

$ {\nabla _{2z}} = {\nabla _{2\eta }} - \left({{\nabla _{2\eta }}z} \right)\left({\frac{{\partial \eta }}{{\partial z}}} \right)\frac{\partial }{{\partial \eta }}, $

(3)

可以得到

$ \frac{\partial }{{\partial t}}\left({{\rho ^ * }\frac{{\partial z}}{{\partial \eta }}} \right) + {\nabla _{2\eta }} \cdot \left({{\rho ^ * }\frac{{\partial z}}{{\partial \eta }}{\mathit{\boldsymbol{V}}_{2z}}} \right) + \frac{\partial }{{\partial \eta }}\left({{\rho ^ * }\frac{{\partial z}}{{\partial \eta }}\dot \eta } \right) = 0, $

(4)

将$\partial z{\rm{/}}\partial \eta = (\partial z{\rm{/}}\partial {p_h})(\partial {p_h}{\rm{/}}\partial \eta) = - \mu {\rm{/}}\rho g$代入公式（4），可得

$ \frac{\partial }{{\partial t}}\left({\frac{{{\rho ^ * }}}{\rho }\mu } \right) + {\nabla _{2\eta }} \cdot \left({\frac{{{\rho ^ * }}}{\rho }\mu {\mathit{\boldsymbol{V}}_{2z}}} \right) + \frac{\partial }{{\partial \eta }}\left({\frac{{{\rho ^ * }}}{\rho }\mu \dot \eta } \right) = 0, $

(5)

定义新变量μ^*，令${\mu ^ * } = ({\rho ^ * }{\rm{/}}\rho)\mu $。将μ^*代入公式（5），可以得到

$ \frac{\partial }{{\partial t}}{({\mu ^ * })_\eta } + {\nabla _{2\eta }} \cdot ({\mu ^ * }{\mathit{\boldsymbol{V}}_{2z}}) + \frac{\partial }{{\partial \eta }}({\mu ^ * }\dot \eta) = 0, $

(6)

公式（6）即为η坐标系中的假不可压连续方程。其中，$\eta = ({p_h} - {p_{ht}}){\rm{/}}\mu $，$\mu = {p_{hs}} - {p_{ht}}$，p_h是气压中的静力分量，p_ht和p_hs分别是模式上、下边界的气压的静力分量。ρ是密度，V_2z是水平风矢量，${\nabla _{2z}} = {\left({\partial /\partial x} \right)_z}\mathit{\boldsymbol{i}} + {\left({\partial /\partial y} \right)_z}\mathit{\boldsymbol{j}}$是z系下的水平梯度算子，${\nabla _{2\eta }} = {\left({\partial /\partial x} \right)_\eta }\mathit{\boldsymbol{i}} + {\left({\partial /\partial y} \right)_\eta }\mathit{\boldsymbol{j}}$是η系下的水平梯度算子。以下，控制方程均投影在η坐标系下，为了简便省略控制方程中部分算子中的角标。

2.2 假不可压控制方程的推导

由Laprise（1992）给出的大气运动控制方程组，并用假不可压连续方程替换其连续方程，可以得到如下形式的假不可压控制方程：

$ \frac{{\partial u}}{{\partial t}} + \left({{\mathit{\boldsymbol{V}}} \cdot \nabla u} \right) + \frac{{{R_d}T}}{p}\frac{{\partial p}}{{\partial x}} + \frac{{\partial p}}{{\partial \eta }}{\left({\frac{{\partial {p_h}}}{{\partial \eta }}} \right)^{ - 1}}\frac{{\partial \phi }}{{\partial x}} = {F_u}, $

(7)

$ \frac{{\partial v}}{{\partial t}} + \left({{\mathit{\boldsymbol{V}}} \cdot \nabla v} \right) + \frac{{{R_d}T}}{p}\frac{{\partial p}}{{\partial y}} + \frac{{\partial p}}{{\partial \eta }}{\left({\frac{{\partial {p_h}}}{{\partial \eta }}} \right)^{ - 1}}\frac{{\partial \phi }}{{\partial y}} = {F_v}, $

(8)

$ \frac{{\partial w}}{{\partial t}} + \left({{\mathit{\boldsymbol{V}}} \cdot \nabla w} \right) - g\left({\frac{1}{\mu }\frac{{\partial p}}{{\partial \eta }} - 1} \right) = {F_w}, $

(9)

$ \frac{{\partial \theta }}{{\partial t}} + \left({{\mathit{\boldsymbol{V}}} \cdot \nabla \theta } \right) = {F_\theta }, $

(10)

$ \frac{\partial }{{\partial t}}{\left({{\mu ^ * }} \right)_\eta } + {\nabla _{2\eta }} \cdot \left({{\mu ^ * }{{\mathit{\boldsymbol{V}}}_{2z}}} \right) + \frac{\partial }{{\partial \eta }}\left({{\mu ^ * }\dot \eta } \right) = 0, $

(11)

$ \frac{{\partial \phi }}{{\partial t}} + \left({{\mathit{\boldsymbol{V}}} \cdot \nabla \phi } \right) - gw = 0, $

(12)

$\mu = \frac{{{\alpha ^ * }}}{\alpha }{\mu ^ * }, $

(13)

${\alpha ^ * } = - \frac{1}{{{\mu ^ * }}}\frac{{\partial \varphi }}{{\partial \eta }}, $

(14)

$ \alpha = - \frac{1}{\mu }\frac{{\partial \varphi }}{{\partial \eta }}, $

(15)

$ p = {p_0}{\left({\frac{{{R_d}\theta }}{{{p_0}\alpha }}} \right)^\gamma }, $

(16)

其中，ϕ = gz是位势，α = 1/ρ是比容，α = 1/ρ^*是假不可压比容，$\gamma = {c_p}{\rm{/}}{c_v} = 1.4$。${\mathit{\boldsymbol{V}}} \cdot \nabla = u\left({\partial {\rm{/}}\partial x} \right) + v\left({\partial /\partial y} \right) + \omega \left({\partial /\partial \eta } \right)$，$\omega = \dot \eta = \left({d\eta {\rm{/}}dt} \right)$。${F_u}, {F_v}, {F_w}, {F_\theta }$是强迫项。

为了与WRF-ARW控制方程的形式一致，将上述假不可压方程进行形式变换，得到了通量形式的控制方程。通量形式的假不可压控制方程形式如下：

$ \frac{{\partial U}}{{\partial t}} + \left({\nabla \cdot {\mathit{\boldsymbol{V}}}u} \right) - \frac{\alpha }{{{\alpha ^ * }}}\frac{\partial }{{\partial x}}\left({p\frac{{\partial \phi }}{{\partial \eta }}} \right) + \frac{\alpha }{{{\alpha ^ * }}}\frac{\partial }{{\partial \eta }}\left({p\frac{{\partial \phi }}{{\partial x}}} \right) = {F_U}, $

(17)

$ \frac{{\partial V}}{{\partial t}} + \left({\nabla \cdot {\mathit{\boldsymbol{V}}}u} \right) - \frac{\alpha }{{{\alpha ^ * }}}\frac{\partial }{{\partial y}}\left({p\frac{{\partial \phi }}{{\partial \eta }}} \right) + \frac{\alpha }{{{\alpha ^ * }}}\frac{\partial }{{\partial \eta }}\left({p\frac{{\partial \phi }}{{\partial y}}} \right) = {F_V}, $

(18)

$ \frac{{\partial W}}{{\partial t}} + \left({\nabla \cdot {\mathit{\boldsymbol{V}}}w} \right) - g\left({\frac{\alpha }{{{\alpha ^ * }}}\frac{{\partial p}}{{\partial \eta }} - {\mu ^ * }} \right) = {F_W}, $

(19)

$ \frac{{\partial \Theta }}{{\partial t}} + \left({\nabla \cdot {\mathit{\boldsymbol{V}}}\theta } \right) = {F_\Theta }, $

(20)

$ \frac{{\partial {\mu ^ * }}}{{\partial t}} + \left({\nabla \cdot {\mathit{\boldsymbol{V}}}} \right) = 0, $

(21)

$ \frac{{\partial \phi }}{{\partial t}} + \frac{1}{{{\mu ^*}}}\left[ {\left({{\mathit{\boldsymbol{V}}} \cdot \nabla \phi } \right) - gW} \right] = 0, $

(22)

$ \mu = \frac{{{\alpha ^ * }}}{\alpha }{\mu ^ * }, $

(23)

$ {\alpha ^ * } = - \frac{1}{{{\mu ^ * }}}\frac{{\partial \varphi }}{{\partial \eta }}, $

(24)

$\alpha = - \frac{1}{\mu }\frac{{\partial \varphi }}{{\partial \eta }}, $

(25)

$p = {p_0}{\left({\frac{{{R_d}\theta }}{{{p_0}\alpha }}} \right)^\gamma }, $

(26)

其中，$U = {\mu ^ * }u$，$V = {\mu ^ * }v$，$W = {\mu ^ * }w$，$\Omega = {\mu ^ * }\dot \eta $，$\Theta = {\mu ^ * }\theta $。$\nabla \cdot \mathit{\boldsymbol{V}}a = \partial Ua/\partial x + \partial Va/\partial y + \partial \Omega a/\partial \eta $，$V \cdot \nabla a = U\left({\partial a{\rm{/}}\partial x} \right) + V\left({\partial a{\rm{/}}\partial y} \right) + \Omega \left({\partial a{\rm{/}}\partial \eta } \right)$，a代表任意气象变量。${F_U}, {F_V}, {F_W}, {F_\Theta }$是强迫项。

2.3 包含水物质的假不可压方程

包含水物质的WRF-ARW控制方程中，连续方程仅描述干空气的质量守恒，同时增加一组水物质预报方程，用于描述水物质的演变（Skamarock et al., 2008）。构建包含水物质的假不可压控制方程时，采用了与WRF模式相同的方法，即连续方程仅针对干空气建立，额外增加一组水物质预报方程。这里略去推导过程，直接给出包含水物质的假不可压方程：

$ \frac{{\partial U}}{{\partial t}} + \left({\nabla \cdot {\mathit{\boldsymbol{V}}}u} \right) + \mu _d^ * \alpha \frac{{\partial p}}{{\partial x}} + \frac{\alpha }{{\alpha _d^ * }}\frac{{\partial p}}{{\partial \eta }}\frac{{\partial \phi }}{{\partial x}} = {F_U}, $

(27)

$ \frac{{\partial V}}{{\partial t}} + \left({\nabla \cdot {\mathit{\boldsymbol{V}}}v} \right) + \mu _d^ * \alpha \frac{{\partial p}}{{\partial y}} + \frac{\alpha }{{\alpha _d^ * }}\frac{{\partial p}}{{\partial \eta }}\frac{{\partial \phi }}{{\partial y}} = {F_V}, $

(28)

$ \frac{{\partial W}}{{\partial t}} + \left({\nabla \cdot {\mathit{\boldsymbol{V}}}w} \right) - g\left({\frac{\alpha }{{\alpha _d^ * }}\frac{{\partial p}}{{\partial \eta }} - \mu _d^ * } \right) = {F_W}, $

(29)

$ \frac{{\partial \Theta }}{{\partial t}} + \left({\nabla \cdot {\mathit{\boldsymbol{V}}}\theta } \right) = {F_\Theta }, $

(30)

$ \frac{{\partial \mu _d^ * }}{{\partial t}} + \left({\nabla \cdot {\mathit{\boldsymbol{V}}}} \right) = 0, $

(31)

$ \frac{{\partial \varphi }}{{\partial t}} + \frac{1}{{\mu _d^ * }}\left[ {\left({{\mathit{\boldsymbol{V}}} \cdot \nabla \varphi } \right) - gW} \right] = 0, $

(32)

$ \frac{{\partial {Q_m}}}{{\partial t}} + \left({\nabla \cdot{\mathit{\boldsymbol{V}}}{q_m}} \right) = {F_{{Q_m}}}, $

(33)

${\mu _d} = \frac{{\alpha _d^ * }}{{{\alpha _d}}}\mu _d^ *, $

(34)

$\alpha _d^ * = - \frac{1}{{\mu _d^ * }}\frac{{\partial \varphi }}{{\partial \eta }}, $

(35)

${\alpha _d} = - \frac{1}{{{\mu _d}}}\frac{{\partial \varphi }}{{\partial \eta }}, $

(36)

$p = {p_0}{\left({\frac{{{R_d}{\theta _m}}}{{{p_0}{\alpha _d}}}} \right)^\gamma }, $

(37)

其中，$\eta = ({p_{dh}} - {p_{dht}}){\rm{/}}{\mu _d}$，${\mu _d} = {p_{dhs}} - {p_{dht}}$，p_dh是干空气中气压的静力分量。$U = \mu _d^ * u$，$V = \mu _d^ * v$，$W = \mu _d^ * w$，$\Omega = \mu _d^ * \dot \eta $，$\Theta = \mu _d^ * \theta $。$\mu _d^ * = ({\alpha _d}{\rm{/}}\alpha _d^ *){\mu _d}$。${\alpha _d} = {\rm{1/}}{\rho _d}$是干空气比容，$\alpha _d^ * = {\rm{1/}}\rho _d^ * $是干空气的假不可压比容，$\rho _d^ * = (\bar \theta {\rm{/}}\theta){\bar \rho _d}$。$\alpha = {\alpha _d}{(1 + {q_v} + {q_c} + {q_r} + {q_i} + ...)^{ - 1}}$，q_v、q_c、q_r、q_i分别是水汽、云滴、雨、冰的混合比。${\alpha ^ * } = \alpha _d^ * {\left({1 + {q_v} + {q_c} + {q_r} + {q_i} + ...} \right)^{ - 1}}$。${\theta _m} = \theta \left[ {1 + \left({{R_v}{\rm{/}}{R_d}} \right){q_v}} \right]$，${Q_m} = \mu _d^ * {q_m}$。${F_{{Q_m}}}$是源汇项。

2.4 扰动形式的假不可压方程

为了减小计算气压梯度项和浮力项的误差，WRF-ARW动力框架使用的是扰动形式的大气运动控制方程。采用相同的方法，将2.3节的假不可压控制方程改写成扰动形式，得到了扰动形式的假不可压控制方程：

$ \begin{array}{l} \frac{{\partial U}}{{\partial t}} + \left({\frac{{\partial Uu}}{{\partial x}} + \frac{{\partial Vu}}{{\partial y}} + \frac{{\partial \Omega u}}{{\partial \eta }}} \right) + \left({\frac{\alpha }{{\alpha _d^ * }}} \right) \cdot \\ \left[ {\mu _d^ * \left({\frac{{\partial \phi '}}{{\partial x}} + \alpha _d^ * \frac{{\partial p'}}{{\partial x}} + \alpha {{_d^ * }^\prime }\frac{{\partial \bar p}}{{\partial x}}} \right) + \frac{{\partial \phi }}{{\partial x}}\left({\frac{{\partial p'}}{{\partial \eta }} - \mu {{_d^ * }^\prime }} \right)} \right] = {F_U}, \end{array} $

(38)

$ \begin{array}{l} \frac{{\partial V}}{{\partial t}} + \left({\frac{{\partial Uv}}{{\partial x}} + \frac{{\partial Vv}}{{\partial y}} + \frac{{\partial \Omega v}}{{\partial \eta }}} \right) + \left({\frac{\alpha }{{\alpha _d^ * }}} \right) \cdot \\ \left[ {\mu _d^ * \left({\frac{{\partial \phi '}}{{\partial y}} + \alpha _d^ * \frac{{\partial p'}}{{\partial y}} + \alpha {{_d^ * }^\prime }\frac{{\partial \bar p}}{{\partial y}}} \right) + \frac{{\partial \phi }}{{\partial y}}\left({\frac{{\partial p'}}{{\partial \eta }} - \mu {{_d^ * }^\prime }} \right)} \right] = {F_V}, \end{array} $

(39)

$ \begin{align} & \frac{\partial W}{\partial t}+\left(\frac{\partial Uw}{\partial x}+\frac{\partial Vw}{\partial y}+\frac{\partial \Omega w}{\partial \eta } \right)-g\frac{\alpha }{\alpha _{d}^{*}}\cdot \\ & \ \ \ \ \ \ \ \ \ \ \ \ \left[ \frac{\partial {p}'}{\partial \eta }-{{{\bar{\mu }}}_{d}}(\frac{\alpha _{d}^{*}}{\alpha }-1) \right]+g\mu {{_{d}^{*}}^{\prime }}={{F}_{W}}, \\ \end{align} $

(40)

$ \frac{{\partial \Theta }}{{\partial t}} + \left({\frac{{\partial U\theta }}{{\partial x}} + \frac{{\partial V\theta }}{{\partial y}} + \frac{{\partial \Omega \theta }}{{\partial \eta }}} \right) = {F_\Theta }, $

(41)

$ \frac{{\partial \mu {{_d^ * }^\prime }}}{{\partial t}} + \left({\frac{{\partial U}}{{\partial x}} + \frac{{\partial V}}{{\partial y}} + \frac{{\partial \Omega }}{{\partial \eta }}} \right) = 0, $

(42)

$ \frac{{\partial \phi '}}{{\partial t}} + \frac{{\rm{1}}}{{\mu _d^ * }}\left[ {\left({U\frac{{\partial \phi }}{{\partial x}} + V\frac{{\partial \phi }}{{\partial y}} + \Omega \frac{{\partial \phi }}{{\partial \eta }}} \right) - gW} \right] = 0, $

(43)

$ \frac{{\partial {Q_m}}}{{\partial t}} + \left({\frac{{\partial U{q_m}}}{{\partial x}} + \frac{{\partial V{q_m}}}{{\partial y}} + \frac{{\partial \Omega {q_m}}}{{\partial \eta }}} \right) = {F_{{Q_m}}}, $

(44)

${\mu '_d} = \frac{{\alpha _d^ * }}{{{\alpha _d}}}\mu _d^ * - {\bar \mu _d}, $

(45)

$ \alpha {{_{d}^{*}}^{\prime }}=-\frac{\rm{1}}{\bar{\mu }_{d}^{*}}\left(\frac{\partial {\phi }'}{\partial \eta }+\alpha _{d}^{*}\mu {{_{d}^{*}}^{\prime }} \right), $

(46)

$ {\alpha '_d} = - \frac{{\rm{1}}}{{{{\bar \mu }_d}}}\left({\frac{{\partial \phi '}}{{\partial \eta }} + {\alpha _d}{{\mu '}_d}} \right), $

(47)

$\frac{{p'}}{{\bar p}} = \gamma \left({\frac{{\Theta '}}{{\bar \Theta }} - \frac{{{{\alpha '}_d}}}{{{{\bar \alpha }_d}}} - \frac{{{{\mu '}_d}}}{{{{\bar \mu }_d}}}} \right), $

(48)

其中，$p = \bar p(\bar z) + p'$，$\phi = \bar \phi (\bar z) + \phi '$，${\alpha _d} = {\bar \alpha _d}(\bar z) + {\alpha '_d}$，$\alpha _{d}^{*}=\bar{\alpha }_{d}^{*}(\bar{z})+\alpha {{_{d}^{*}}^{\prime }}$，${\mu _d} = {\bar \mu _d}(x, y) + {\mu '_d}$，$\mu _{d}^{*}=\bar{\mu }_{d}^{*}(x, y)+\mu {{_{d}^{*}}^{\prime }}$。

2.5 时间离散的假不可压方程

WRF模式采用时步分离的积分方案，即用较大的时间步长积分慢波，用时步较小的声模循环积分高频波，兼顾计算效率和稳定。假不可压模式中，同样采用了时步分离的积分方案。这里直接给出时间离散形式的假不可压方程：

$ \begin{array}{l} \frac{{\delta U''}}{{\delta \tau }} + \left({\frac{{{\alpha ^{{t^*}}}}}{{\alpha _d^{ * {t^*}}}}} \right)\left[ {\mu {{_d^*}^{{t^ * }}}\left({\alpha {{_d^*}^{{t^ * }}}\frac{{\partial {{p''}^\tau }}}{{\partial x}} + \alpha {{_d^*}^{\prime \prime }}^\tau \frac{{\partial \bar p}}{{\partial x}} + \frac{{\partial {{\phi ''}^\tau }}}{{\partial x}}} \right) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\frac{{\partial {\phi ^{{t^*}}}}}{{\partial x}}{{\left({\frac{{\partial p''}}{{\partial \eta }} - \mu {{_d^*}^{\prime \prime }}} \right)}^\tau }} \right] = R_U^{{t^*}}, \end{array} $

(49)

$ \begin{array}{l} \frac{{\delta V''}}{{\delta \tau }} + \left({\frac{{{\alpha ^{{t^*}}}}}{{\alpha _d^{ * {t^*}}}}} \right)\left[ {\mu {{_d^*}^{{t^ * }}}\left({\alpha {{_d^*}^{{t^ * }}}\frac{{\partial {{p''}^\tau }}}{{\partial y}} + \alpha {{_d^*}^{\prime \prime }}^\tau \frac{{\partial \bar p}}{{\partial y}} + \frac{{\partial {{\phi ''}^\tau }}}{{\partial y}}} \right) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{ }}\left. {\frac{{\partial {\phi ^{{t^*}}}}}{{\partial y}}{{\left({\frac{{\partial p''}}{{\partial \eta }} - \mu {{_d^*}^{\prime \prime }}} \right)}^\tau }} \right] = R_V^{{t^*}}, \end{array} $

(50)

$ \frac{{\delta \mu {{_d^ * }^{\prime \prime }}}}{{\delta \tau }} + {\left({\frac{{\partial {U^{\prime \prime }}}}{{\partial x}} + \frac{{\partial {V^{\prime \prime }}}}{{\partial y}} + \frac{{\partial {\Omega ^{\prime \prime }}}}{{\partial \eta }}} \right)^{\tau + \Delta \tau }} = R_{\mu _d^ * }^{{t^*}}, $

(51)

$ \frac{{\delta \Theta ''}}{{\delta \tau }} + {\left({\frac{{\partial U''{\theta ^{{t^*}}}}}{{\partial x}} + \frac{{\partial V''{\theta ^{{t^*}}}}}{{\partial y}} + \frac{{\partial \Omega ''{\theta ^{{t^*}}}}}{{\partial \eta }}} \right)^{\tau + \Delta \tau }} = R_\Theta ^{{t^*}}, $

(52)

$ \begin{array}{l} \frac{{\delta W''}}{{\delta \tau }} - \\ g{\overline {\left\{ {{{\left({\frac{\alpha }{{\alpha _d^ * }}} \right)}^{{t^ * }}}\left[ {\frac{\partial }{{\partial \eta }}\left({C\frac{{\partial \phi ''}}{{\partial \eta }}} \right) + \frac{\partial }{{\partial \eta }}\left({\frac{{c_s^2}}{{\alpha _d^{{t^*}}}}\frac{{\Theta ''}}{{{\Theta ^{{t^*}}}}}} \right)} \right] - \mu {{_d^*}^{\prime \prime }}} \right\}} ^\tau } = R_W^{{t^ * }}, \end{array} $

(53)

$ \frac{\delta {\phi }''}{\delta \tau }+\frac{1}{\mu {{_{d}^{*}}^{{{t}^{*}}}}}\left({{{{\Omega }''}}^{\tau +\Delta \tau }}\frac{\partial {{\phi }^{{{t}^{*}}}}}{\partial \eta }-g{{{\bar{{W}''}}}^{\tau }} \right)=R_{\phi }^{{{t}^{*}}}, $

(54)

$ {{{\mu }''}_{d}}={{\left(\frac{\alpha _{d}^{*}}{{{\alpha }_{d}}} \right)}^{{{t}^{*}}}}\mu {{_{d}^{*}}^{\prime\prime }}, $

(55)

$ \alpha {{_{d}^{*}}^{\prime\prime }}=-\frac{1}{\mu {{_{d}^{*}}^{{{t}^{*}}}}}\left(\frac{\partial {\phi }''}{\partial \eta }+\alpha {{_{d}^{*}}^{{{t}^{*}}}}\mu {{_{d}^{*}}^{\prime\prime }} \right), $

(56)

$ {{{\alpha }''}_{d}}=-\frac{1}{\mu _{d}^{{{t}^{*}}}}\left(\frac{\partial {\phi }''}{\partial \eta }+\alpha _{d}^{{{t}^{*}}}{{{{\mu }''}}_{d}} \right), $

(57)

$ \frac{{{p}''}}{{{p}^{{{t}^{*}}}}}=\gamma \left(\frac{{{\Theta }''}}{{{\Theta }^{{{t}^{*}}}}}-\frac{{{{{\alpha }''}}_{d}}}{\alpha _{d}^{{{t}^{*}}}}-\frac{{{{{\mu }''}}_{d}}}{\mu _{d}^{{{t}^{*}}}} \right), $

(58)

方程右端项的离散形式如下：

$ \begin{align} & R_{U}^{{{t}^{*}}}=-\left(\frac{\partial Uu}{\partial x}+\frac{\partial Vu}{\partial y}+\frac{\partial \Omega u}{\partial \eta } \right)-\left(\frac{\alpha }{\alpha _{d}^{*}} \right)\cdot \\ & \ \ \left[ \mu _{d}^{*}\left(\frac{\partial {\phi }'}{\partial x}+\alpha _{d}^{*}\frac{\partial {p}'}{\partial x}+\alpha {{_{d}^{*}}^{\prime }}\frac{\partial \bar{p}}{\partial x} \right)+\frac{\partial \phi }{\partial x}\left(\frac{\partial {p}'}{\partial \eta }-\mu {{_{d}^{*}}^{\prime }} \right) \right]+{{F}_{U}}, \\ \end{align} $

(59)

$ \begin{align} & R_{V}^{{{t}^{*}}}=-\left(\frac{\partial Uv}{\partial x}+\frac{\partial Vv}{\partial y}+\frac{\partial \Omega v}{\partial \eta } \right)-\left(\frac{\alpha }{\alpha _{d}^{*}} \right)\cdot \\ & \ \left[ \mu _{d}^{*}\left(\frac{\partial {\phi }'}{\partial y}+\alpha _{d}^{*}\frac{\partial {p}'}{\partial y}+\alpha {{_{d}^{*}}^{\prime }}\frac{\partial \bar{p}}{\partial y} \right)+\frac{\partial \phi }{\partial y}\left(\frac{\partial {p}'}{\partial \eta }-\mu {{_{d}^{*}}^{\prime }} \right) \right]+{{F}_{V}}, \\ \end{align} $

(60)

$ \begin{align} & R_{W}^{{{t}^{*}}}=-\left(\frac{\partial Uw}{\partial x}+\frac{\partial Vw}{\partial y}+\frac{\partial \Omega w}{\partial \eta } \right)+g\left(\frac{\alpha }{\alpha _{d}^{*}} \right)\cdot \\ & \ \ \ \ \left[ \frac{\partial {p}'}{\partial \eta }-{{{\bar{\mu }}}_{d}}\left(\frac{\alpha _{d}^{*}}{\alpha }-1 \right) \right]-\mu {{_{d}^{*}}^{\prime }}g+{{F}_{W}}, \\ \end{align} $

(61)

$ R_{\mu _{d}^{*}}^{{{t}^{*}}}=-\left(\frac{\partial U}{\partial x}+\frac{\partial V}{\partial y}+\frac{\partial \Omega }{\partial \eta } \right), $

(62)

$ R_{\Theta }^{{{t}^{*}}}=-\left(\frac{\partial U\theta }{\partial x}+\frac{\partial V\theta }{\partial y}+\frac{\partial \Omega \theta }{\partial \eta } \right)+{{F}_{\Theta }}, $

(63)

$ R_{\phi }^{{{t}^{*}}}=-\frac{1}{\mu _{d}^{*}}\left(U\frac{\partial \phi }{\partial x}+V\frac{\partial \phi }{\partial y}+\Omega \frac{\partial \phi }{\partial \eta }-gW \right), $

(64)

其中，${U}''=U-{{U}^{{{t}^{*}}}}$，${V}''=V-{{V}^{{{t}^{*}}}}$，${W}''=W-{{W}^{{{t}^{*}}}}$，${\Omega }''=\Omega -{{\Omega }^{{{t}^{*}}}}$，${\Theta }''=\Theta -{{\Theta }^{{{t}^{*}}}}$，${\varphi }''={\varphi }'-{{{\varphi }'}^{{{t}^{*}}}}$，${{{\alpha }''}_{d}}={{{\alpha }'}_{d}}-{{{\alpha }'}_{d}}^{{{t}^{*}}}$，$\alpha {{_{d}^{*}}^{\prime\prime }}=\alpha {{_{d}^{*}}^{\prime }}-\alpha {{_{d}^{*}}^{\prime }}^{{{t}^{*}}}$，${{{\mu }''}_{d}}={{{\mu }'}_{d}}-{{{\mu }'}_{d}}^{{{t}^{*}}}$，$\mu {{_{d}^{*}}^{\prime\prime }}=\mu {{_{d}^{*}}^{\prime }}-\mu {{_{d}^{*}}^{\prime }}^{{{t}^{*}}}$，$c_{s}^{2}=\gamma {{p}^{{{t}^{*}}}}\alpha _{d}^{{{t}^{*}}}$，$C=c_{s}^{2}\rm{/(}\mu _{d}^{{{t}^{*}}}\alpha {{_{d}^{{{t}^{*}}}}^{2}}\rm{)}$。

时间离散算子为

$ \frac{\delta a}{\delta \tau }=\frac{{{a}^{\tau +\Delta \tau }}-{{a}^{\tau }}}{\Delta \tau }, $

(65)

${{\bar{a}}^{\tau }}=\frac{1+\beta }{2}{{a}^{\tau +\Delta \tau }}+\frac{1-\beta }{2}{{a}^{\tau }}, $

(66)

其中，Δτ为声模循环中的时间步长，β由用户给定。

2.6 空间离散的假不可压方程

将2.5节得到的假不可压方程投影到Arakawa C交错网格上，即可得到空间离散形式的假不可压方程：

$ \begin{align} & \frac{\delta {U}''}{\delta \tau }+{{\overline{\left(\frac{{{\alpha }^{{{t}^{*}}}}}{\alpha {{_{d}^{*}}^{{{t}^{*}}}}} \right)}}^{x}}\left[ {{\overline{\mu {{_{d}^{*}}^{{{t}^{*}}}}}}^{x}}\left({{\overline{\alpha {{_{d}^{*}}^{{{t}^{*}}}}}}^{x}}\frac{\delta {{{{p}''}}^{\tau }}}{\delta x}+{{\overline{\alpha {{_{d}^{*}}^{\prime\prime }}^{\tau }}}^{x}}\frac{\delta \bar{p}}{\delta x}+\frac{\delta {{\overline{{{{{\phi }''}}^{\tau }}}}^{\eta }}}{\delta x} \right)+ \right. \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left. \frac{\delta {{\overline{{{\phi }^{{{t}^{*}}}}}}^{\eta }}}{\delta x}{{\left(\frac{\delta {{\overline{{{\overline{{{p}''}}}^{x}}}}^{\eta }}}{\delta \eta }-{{\overline{\mu {{_{d}^{*}}^{\prime\prime }}}}^{x}} \right)}^{\tau }} \right]=R_{U}^{{{t}^{*}}}, \\ \end{align} $

(67)

$ \begin{align} & \frac{\delta {V}''}{\delta \tau }+{{\overline{\left(\frac{{{\alpha }^{{{t}^{*}}}}}{\alpha {{_{d}^{*}}^{{{t}^{*}}}}} \right)}}^{y}}\left[ {{\overline{\mu {{_{d}^{*}}^{{{t}^{*}}}}}}^{y}}\left({{\overline{\alpha {{_{d}^{*}}^{{{t}^{*}}}}}}^{y}}\frac{\delta {{{{p}''}}^{\tau }}}{\delta y}+{{\overline{\alpha {{_{d}^{*}}^{\prime\prime }}^{\tau }}}^{y}}\frac{\delta \bar{p}}{\delta y}+\frac{\delta {{\overline{{{{{\phi }''}}^{\tau }}}}^{\eta }}}{\delta y} \right)+ \right. \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left. \frac{\delta {{\overline{{{\phi }^{{{t}^{*}}}}}}^{\eta }}}{\delta y}{{\left(\frac{\delta {{\overline{{{\overline{{{p}''}}}^{y}}}}^{\eta }}}{\delta \eta }-{{\overline{\mu {{_{d}^{*}}^{\prime\prime }}}}^{y}} \right)}^{\tau }} \right]=R_{V}^{{{t}^{*}}}, \\ \end{align} $

(68)

$ \begin{align} & \frac{\delta {W}''}{\delta \tau }-g\cdot \\ & {{\overline{\left\{ {{\overline{{{\left(\frac{\alpha }{\alpha _{d}^{*}} \right)}^{{{t}^{*}}}}}}^{\eta }}\left[ \frac{\delta }{\delta \eta }\left(C\frac{\delta {\phi }''}{\delta \eta } \right)+\frac{\delta }{\delta \eta }\left(\frac{c_{s}^{2}}{\alpha _{d}^{{{t}^{*}}}}\frac{{{\Theta }''}}{{{\Theta }^{{{t}^{*}}}}} \right) \right]-\mu {{_{d}^{*}}^{\prime\prime }} \right\}}}^{\tau }}=R_{W}^{{{t}^{*}}}, \\ \end{align} $

(69)

$ \frac{\delta \mu {{_{d}^{*}}^{\prime\prime }}}{\delta \tau }+{{\left(\frac{\delta {U}''}{\delta x}+\frac{\delta {V}''}{\delta y}+\frac{\delta {\Omega }''}{\delta \eta } \right)}^{\tau +\Delta \tau }}=R_{\mu _{d}^{*}}^{{{t}^{*}}}, $

(70)

$ \frac{\delta {\Theta }''}{\delta \tau }+{{\left(\frac{\delta {U}''{{\overline{{{\theta }^{{{t}^{*}}}}}}^{x}}}{\delta x}+\frac{\delta {V}''{{\overline{{{\theta }^{{{t}^{*}}}}}}^{y}}}{\delta y}+\frac{\delta {\Omega }''{{\overline{{{\theta }^{{{t}^{*}}}}}}^{\eta }}}{\delta \eta } \right)}^{\tau +\Delta \tau }}=R_{\Theta }^{{{t}^{*}}}, $

(71)

$ \frac{\delta {\phi }''}{\delta \tau }+\frac{1}{{{\overline{\mu {{_{d}^{*}}^{{{t}^{*}}}}}}^{\eta }}}\left({{{{\Omega }''}}^{\tau +\Delta \tau }}\frac{\delta {{\overline{{{\phi }^{{{t}^{*}}}}}}^{\eta }}}{\delta \eta }-g{{{\bar{{W}''}}}^{\tau }} \right)=R_{\phi }^{{{t}^{*}}}, $

(72)

$ {{{\mu }''}_{d}}={{\left(\frac{\alpha _{d}^{*}}{{{\alpha }_{d}}} \right)}^{{{t}^{*}}}}\mu {{_{d}^{*}}^{\prime\prime }}, $

(73)

$ \alpha {{_{d}^{*}}^{\prime\prime }}=-\frac{1}{\mu {{_{d}^{*}}^{{{t}^{*}}}}}\left(\frac{\delta {\phi }''}{\delta \eta }+\alpha {{_{d}^{*}}^{{{t}^{*}}}}\mu {{_{d}^{*}}^{\prime\prime }} \right), $

(74)

$ {{{\alpha }''}_{d}}=-\frac{1}{\mu _{d}^{{{t}^{*}}}}\left(\frac{\delta {\phi }''}{\delta \eta }+\alpha _{d}^{{{t}^{*}}}{{{{\mu }''}}_{d}} \right), $

(75)

$\frac{{{p}''}}{{{p}^{{{t}^{*}}}}}=\gamma \left(\frac{{{\Theta }''}}{{{\Theta }^{{{t}^{*}}}}}-\frac{{{{{\alpha }''}}_{d}}}{\alpha _{d}^{{{t}^{*}}}}-\frac{{{{{\mu }''}}_{d}}}{\mu _{d}^{{{t}^{*}}}} \right), $

(76)

方程右端项的离散形式如下：

$ \begin{align} & R_{U}^{{{t}^{*}}}=-{{\overline{\left(\frac{\alpha }{\alpha _{d}^{*}} \right)}}^{x}}\left[ {{\overline{\mu _{d}^{*}}}^{x}}\left({{\overline{\alpha _{d}^{*}}}^{x}}\frac{\delta {p}'}{\delta x}+{{\overline{\alpha {{_{d}^{*}}^{\prime }}}}^{x}}\frac{\delta \bar{p}}{\delta x}+\frac{\delta {{\overline{{{\phi }'}}}^{\eta }}}{\delta x} \right)+ \right. \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left. \frac{\delta {{\overline{\phi }}^{\eta }}}{\delta x}\left(\frac{\delta {{\overline{{{\overline{{{p}'}}}^{x}}}}^{\eta }}}{\delta \eta }-{{\overline{\mu {{_{d}^{*}}^{\prime }}}}^{x}} \right) \right]+{{F}_{U}}, \\ \end{align} $

(77)

$ \begin{align} & R_{V}^{{{t}^{*}}}=-{{\overline{\left(\frac{\alpha }{\alpha _{d}^{*}} \right)}}^{y}}\left[ {{\overline{\mu _{d}^{*}}}^{y}}\left({{\overline{\alpha _{d}^{*}}}^{y}}\frac{\delta {p}'}{\delta y}+{{\overline{\alpha {{_{d}^{*}}^{\prime }}}}^{y}}\frac{\delta \bar{p}}{\delta y}+\frac{\delta {{\overline{{{\phi }'}}}^{\eta }}}{\delta y} \right)+ \right. \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left. \frac{\delta {{\overline{\phi }}^{\eta }}}{\delta y}\left(\frac{\delta {{\overline{{{\overline{{{p}'}}}^{y}}}}^{\eta }}}{\delta \eta }-{{\overline{\mu {{_{d}^{*}}^{\prime }}}}^{y}} \right) \right]+{{F}_{V}}, \\ \end{align} $

(78)

$ R_{W}^{{{t}^{*}}}=g{{\overline{\left(\frac{\alpha }{\alpha _{d}^{*}} \right)}}^{\eta }}\left[ \frac{\delta {p}'}{\delta \eta }-\bar{\mu }_{d}^{*}\left(\frac{\alpha _{d}^{*}}{\alpha }-1 \right) \right]-\mu {{_{d}^{*}}^{\prime }}g+{{F}_{W}}, $

(79)

$ R_{\mu _{d}^{*}}^{{{t}^{*}}}=-\left(\frac{\delta U}{\delta x}+\frac{\delta V}{\delta y}+\frac{\delta \Omega }{\delta \eta } \right), $

(80)

$ R_{\Theta }^{{{t}^{*}}}=-\left(\frac{\delta U\theta }{\delta x}+\frac{\delta V\theta }{\delta y}+\frac{\delta \Omega \theta }{\delta \eta } \right)+{{F}_{\Theta }}, $

(81)

$ R_{\phi }^{{{t}^{*}}}=-\frac{1}{{{\overline{\mu _{d}^{*}}}^{\eta }}}\left(U\frac{\delta \phi }{\delta x}+V\frac{\delta \phi }{\delta y}+\Omega \frac{\delta \phi }{\delta \eta }-gW \right), $

(82)

空间离散算子如下：

$ \frac{\delta a}{\delta x}=\frac{{{a}_{i+1/2}}-{{a}_{i-1/2}}}{\Delta x}, $

(83)

$ \frac{\delta a}{\delta y}=\frac{{{a}_{j+1/2}}-{{a}_{j-1/2}}}{\Delta y}, $

(84)

$ \frac{{\delta a}}{{\delta \eta }} = \frac{{{a_{k + 1/2}} - {a_{k - 1/2}}}}{{\Delta \eta }}, $

(85)

${\bar a^\eta }{|_{k + 1/2}} = \frac{1}{2}\left({\frac{{\Delta {\eta _k}}}{{\Delta {\eta _{k + 1/2}}}}{a_{k + 1}} + \frac{{\Delta {\eta _{k + 1}}}}{{\Delta {\eta _{k + 1/2}}}}{a_k}} \right), $

(86)

其中，i, j, k分别是x, y, _η轴的格点坐标。

3 数值试验 3.1 湿热泡对流试验 3.1.1 数值试验的设计方案

利用WRF模式大涡模拟模块（WRF-LES）和假不可压模式大涡模拟模块（Pseudo-LES）分别进行湿热泡对流试验。模拟区域的设置：水平0≤x≤7800 m，0≤y≤7800 m，格距200 m；垂直方向h≤z≤10000 m；垂直方向分为50层。时间步长为1 s，模拟时间长度为30 min。微物理方案选用Lin方案（Lin et al., 1983），该方案包含了水汽、云水、雨、云冰、雪、霰六种水物质，较为复杂精细，适用于高分辨率模拟。探空廓线采用的是Yau（1980）的湿廓线，该廓线水平风速为零，低层相对湿度大，有利于对流的发展。在初始场设置了一个三维热泡，作为对流的启动机制（参考Klein, 2009）。热泡的半径（r）及扰动位温（θ）的给定如下：

$\theta = \left\{ \begin{array}{l} 3{\cos ^2}\left({\frac{\pi }{2}r} \right){\rm{ }}(r \le 1), \\ 0, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他, \end{array} \right.$

(87)

其中，

$r = \sqrt {{{\left({\frac{{x - 3800}}{{2000}}} \right)}^2} + {{\left({\frac{{y - 3800}}{{2000}}} \right)}^2} + {{\left({\frac{{z - 1200}}{{600}}} \right)}^2}, } $

(88)

上式中，x、y、z代表三个坐标轴，距离的单位是m。

假不可压模式建立在η坐标系下，而η坐标系具备处理复杂地形的能力。为了检验假不可压模式对非平坦地形个例的模拟效果，探究不同地形条件下湿热泡发展演变的差异，并考察假不可压模式与原模式模拟结果的相似性，我们在模拟区域内分别设置了平缓山丘和陡峭山丘，进行对照数值试验，比较不同地形下假不可压模式的模拟结果，并对比假不可压模式与WRF模式的模拟结果。两种地形条件的山丘高度（h）的给定如下，

平缓山丘：

$ h\left({x, y} \right) = \left\{ \begin{array}{l} 200\cos \left[ {\frac{\pi }{{1200}}\left({x - 3800} \right)} \right]\\ \cos \left[ {\frac{\pi }{{1200}}\left({y - 3800} \right)} \right], {\rm{ }}\begin{array}{*{20}{c}} {3200 \le x \le 4400}\\ {3200 \le y \le 4400} \end{array}, \\ 0, \;\;\;\;\;\;\;其他, \end{array} \right. $

(89)

陡峭山丘：

$h\left({x, y} \right) = \left\{ \begin{array}{l} 400\cos \left[ {\frac{{\rm{ \mathsf{ π} }}}{{1200}}\left({x - 3800} \right)} \right]\\ \cos \left[ {\frac{{\rm{ \mathsf{ π} }}}{{1200}}\left({y - 3800} \right)} \right], {\rm{ }}\begin{array}{*{20}{c}} {3200 \le x \le 4400}\\ {3200 \le y \le 4400} \end{array}, \\ 0, \;\;\;\;\;\;\;其他, \end{array} \right.$

(90)

其中，h是地形高度，距离的单位是m。

3.1.2 模拟结果及分析

（a）平缓山丘条件下模拟结果及分析

初始场中的湿热泡在浮力的驱动下上升，上升过程中膨胀冷却降温，饱和度增加，达饱和后凝结成云。图 1、图 2分别是第11 min、第13 min云水含量的垂直剖面，第11 min时对流云处于发展阶段，云底高度约为1.2 km，云顶伸展到接近5 km，云水含量最大值超过了3 g kg⁻¹。2 min后积云发展到强盛阶段，云顶伸展至6 km，云水含量明显减少，最大值小于2 g kg^-1。图 3是云水含量的高度—时间分布，两种模式对云体的模拟结果一致，仅在细节刻画上有细微区别。图 4、图 5分别是第11 min、第13 min垂直速度的分布，由图可见对流发展十分旺盛，上升速度的最大值超过了20 m s^-1。进一步分析微物理场和动力场（图 1、图 4）可以看出，对流发展阶段云内上升气流很强，这有利于小云滴迅速凝结增长，云水含量增加；而云滴凝结释放潜热，增大浮力，又会进一步促进对流发展，形成正反馈效应；同时，在垂直对流的作用下大云滴会集中于云的上部，因此云水含量大值区出现在云的上部。图 6是垂直速度的高度—时间分布，如图所示对流热泡经历了两次发展，第一次发展是从初始时刻至第17 min，第二次发展则是在模拟的第17至24 min，两种模式均模拟出了对流热泡的两次发展过程。总体而言，假不可压模式准确的还原了WRF模式的动力场。由于假不可压模式没有对动量方程做近似，完整地保留了驱动对流云发展的浮力项与扰动气压梯度力项（Sun et al., 2012），因此假不可压模式的模拟结果理应与WRF模式一致。

（b）陡峭山丘条件下模拟结果及分析

云水含量的垂直剖面，与平缓地形相比，增高地形后相同时刻对流云发展的更高，第11 min时云顶高度已经超过5 km，13 min时达到6.5 km，云水含量与平缓地形下的模拟结果大致相同。图 9给出了增高地形后云水含量的高度—时间分布，可以看出在陡峭山丘条件下，假不可压模式的模拟结果仍与原模式的模拟结果非常接近。图 10、图 11分别是第11 min、第13 min垂直速度分布，与平缓地形相比，抬高地形后对流发展的更为强烈，第13 min时对流云内最大垂直速度超过了24 m s^-1。比较假不可压模式和WRF模式模拟出的垂直速度随时间的演变（图 12），可以看出在陡峭山丘条件下假不可压模式模拟的动力场与WRF模式的模拟结果基本一致。

综合分析两种地形条件下的模拟结果，可以看出建立在η坐标系下的假不可压模式能够模拟出非平坦地形条件下湿热泡的发展演变，且两种地形条件下假不可压模式的模拟结果均接近WRF模式的模拟结果。由于假不可压模式未对动量方程做任何近似，并且与WRF模式采用相同的微物理方案，因此这一结果是合理可信的。另外，在本次数值试验中湿热泡的发展演变对于地形的变化是敏感的，增加地形高度加速了对流云的发展。由于湿热泡理想试验不考虑下垫面的热力效应，因此是地形的动力效应导致对流发展速度的变化。

3.2 重力流数值试验 3.2.1 数值试验的设计方案

利用WRF模式二维重力流模块（WRF-grav2d_ x）和假不可压模式二维重力流模块（Pseudo-grav2d_x）分别进行二维重力流数值试验（参考Straka et al., 1993）。模拟区域的设置是x方向0≤x≤40000 m，格距100 m；垂直方向0≤z≤6409 m，分为65层。时间步长为1 s，模拟时间长度为15 min。探空廓线使用的是WRF模式提供的静风廓线。在初始场3000 m高度处设置了一个二维的冷泡，冷泡的半径（r）及扰动位温（θ）的给定如下：

$\theta = \left\{ \begin{array}{l} - 15{\left({\frac{{1000}}{p}} \right)^{\frac{R}{{{c_p}}}}}\frac{{\cos (\pi r) + 1}}{{\rm{2}}}{\rm{ }}(r \le 1), \\ 0, \;\;\;\;\;\;\;其他, \end{array} \right.$

(91)

其中，

$r = \sqrt {{{\left({\frac{{x - 19900}}{{4000}}} \right)}^2} + {{\left({\frac{{z - 3000}}{{2000}}} \right)}^2}, } $

(92)

上式中，x是水平坐标轴，z是垂直坐标轴，距离的单位是m。

3.2.2 模拟结果及分析

图 13是初始时刻扰动位温的垂直剖面，如图所示在3 km的高度设置了一个椭圆形的冷泡，长轴长8 km，短轴长4 km，同一高度冷泡中心的扰动位温低于两侧。冷泡受到的浮力小于重力，在重力的作用下冷泡将下沉，同时冷泡的形状也会扭曲。冷泡重力流试验能够检验假不可压模式对于斜压非线性流体的模拟能力。图 14给出了第5 min扰动位温的空间分布，此时冷泡已经落到地面，在地面的约束作用下，冷空气向两侧推进，推开周围的暖空气。随后冷空气顶部边界处会形成Kelvin-Helmholz切变不稳定涡旋（Straka et al., 1993）。第10 min时（图 15），切变涡旋已经具有闭合结构，冷涡旋的外侧有一个冷空气塔。接下来冷空气会继续向两侧推进，涡旋范围扩大但是强度减弱。图 16是第15 min扰动位温的垂直剖面，由图可见此时涡旋结构仍十分清晰，但中心位温较最初已经明显升高。图 17是扰动位温的高度—时间分布，与WRF模式模拟结果相比可知假不可压模式准确的模拟出了冷泡的演变。二维重力流数值试验证明了假不可压模式具备准确模拟斜压非线性流体发展演变的能力。

4 总结与讨论 4.1 总结

建立声波滤除模式可以消除声波对积分时步的严格限制，从而可以用较大的时间步长进行长时间的稳定积分，节约了计算资源，提高了计算效率。经典的滤波控制方程有包辛尼斯克方程、滞弹性方程和假不可压方程。现有的研究表明包辛尼斯克方程不适用于研究深对流；滞弹性方程则需要假设位温的层结分布，适用范围有局限性；假不可压方程不仅更接近真实大气，而且成立条件容易满足，适合作为控制方程，用于构建数值预报模式。因此本文采用了假不可压近似理论作为理论基础，进行滤除声波模式的构建。

WRF模式广泛用于数值预报和研究，拥有成熟的平流方案和微物理模块。WRF-ARW动力框架建立在地形追随坐标系下，可以处理复杂地形问题。本文选取WRF-ARW动力框架作为模式的框架基础，用推导出的假不可压控制方程改写WRF模式的相关代码，成功构建了假不可压模式。建立假不可压模式的步骤如下：

（1）将假不可压连续方程投影到η坐标系，并改写成与WRF-ARW动力框架中的连续方程相近的形式。在此基础上，根据WRF-ARW控制方程的形式，推导出了一组通量形式的假不可压大气运动控制方程，作为假不可压模式的控制方程。

（2）推导出包含水物质的假不可压控制方程，将控制方程改写成扰动形式。对控制方程进行时间和空间离散，得到离散形式的假不可压控制方程。

（3）采用与WRF-ARW相近的数值求解方案：时间积分方案沿用WRF-ARW的原方案；控制方程的求解顺序和求解方法也与WRF-ARW相近。

（4）依据假不可压理论，改写WRF-ARW动力框架部分的代码，改写大涡模拟模块和二维重力流模块的初始化部分的代码，得到了假不可压大涡模式和假不可压二维重力流模式。

为了检验假不可压模式的准确性，本文进行了两组对照数值试验。首先，用Pseudo-LES和WRF-LES分别进行了湿热泡对流试验，试验的设计方案一致。分析第一组对照试验的模拟结果，可知假不可压模式模拟的湿热泡的演变与WRF模式模拟的结果非常接近；假不可压模式具备模拟非平坦地形的个例的能力，且对地形的变化是敏感的。其次，用Pseudo-grav2d_x和WRF-grav2d_x采用相同的设计方案分别进行重力流试验。分析第二组对照试验的模拟结果，可以看出假不可压模式模拟出的冷泡的演变和WRF模式的模拟结果几乎完全一致，表明假不可压模式能够准确模拟斜压非线性流体在重力作用下的演变过程。两组对照试验中，假不可压模式均准确还原了WRF模式的模拟结果，这说明在WRF模式的模式框架下建立的假不可压模式是合理可信的。

4.2 讨论

理论上讲，本文构建的假不可压模式的控制方程中已经滤除了声波，因此在数值求解时也就不需要使用声模循环来积分高频波，并且模式对时间步长的限制也没有全弹性模式苛刻。然而，在实践中我们发现与纯理论的滤波方程相比，在数值模式中滤除全部的高频波是一个更加复杂的系统性工程，仅采用滤波方程不足以彻底滤除数值解中的高频波。若想彻底滤除数值解中的高频波，实现增大模拟的时间步长提高计算效率的目标，还有如下问题需要解决：WRF模式滤波方案中包括了滤除声波和外部模态的算法，而假不可压模式仅需滤除外部模态，因此需要设计新的算法滤除假不可压模式数值解中的外部模态；目前建立的假不可压模式的数值求解方案包含了声模循环，去除声模循环是提高计算效率的关键。

在WRF-ARW模式框架下引入假不可压方程，建立了假不可压模式，并通过对照数值试验验证了假不可压模式的准确性，这是本文的主要创新点；尚未从根本上改变WRF-ARW数值求解方案的架构，并设计出新的滤波算法，导致构建的假不可压模式无法使用更大的时间步长进行数值模拟，体现假不可压模式的节约计算资源优势，这是当前工作的不足。在后面的研究中，将着力于改进假不可压模式的数值求解方案和滤波算法，去除声模循环，尝试增大模式的时间步长，最终实现提高模式的模拟速度、节约计算资源的目标。