1 引言

数值天气预报是业务天气预报的重要技术支撑，而大气数值模式的发展是数值天气预报的基础。尽可能提高模式水平分辨率是提高数值天气预报准确率的一个重要方面。目前许多数值天气预报中心的业务模式水平分辨率已达到2~10 km（Davies et al., 2005; Saito et al., 2006; Charles and Colle, 2009; 陈敏等，2011；Baldauf et al., 2011; Seity et al., 2011；陈葆德等，2013；Benjamin et al., 2016）。Hong and Dudhia（2012）预计在未来几年内一些数值天气预报中心的模式水平分辨率可达1 km，甚至1 km以下。模式向高分辨率发展不仅在计算上对模式动力框架提出了更高的要求，而且需要在物理上对各种参数化方案进行重新思考与认识。其中以湍流与边界层过程为代表的高分辨率下次网格尺度湍流混合的参数化问题是最重要的方面之一（Dudhia, 2014; 刘辉志等，2013）。

从尺度与行星边界层（Planetary Boundary Layer，简称PBL）高度相当的大涡，到尺度为毫米量级的可耗散湍流能量的涡动，大气湍流具有连续的、不同的空间尺度。而在数值模式中，需要给定网格距，因此必须对连续的尺度进行截断。大气运动中含能湍涡的长度尺度l_e在深湿对流和对流边界层内可达公里量级（Shin and Hong, 2015）。

当模式网格距Δ远大于含能湍涡的长度尺度l_e时，例如，传统的分辨率较粗的中尺度数值预报模式（水平分辨率10~100 km），通常可认为大气边界层内所有湍流运动均为次网格尺度，需要对所有湍流进行参数化。此时网格内具有足够多的湍流统计样本，参数化需要表征的是湍流样本集合平均的效应，网格空间内具有相同的统计特征，满足水平均一性条件，所以主要考虑垂直方向上的湍流混合作用。因此，在中尺度模式中，通常使用一维大气边界层参数化方案（以下简称PBL方案）来表征次网格湍流混合的作用，如局地MYJ（Mellor-Yamada-Janjic）方案（Janjić, 2001）和非局地YSU（Yonsei University）方案（Hong et al., 2006）。

当模式网格距Δ远小于含能湍涡的长度尺度l_e时，经常使用大涡模式（Large-Eddy Simulation，简称LES）进行湍流模拟。LES的水平网格距通常要求小于100 m，以确保分辨率落在惯性子区，此时模式动力过程能够直接解析尺度较大的含能湍涡，仅对小尺度涡动进行参数化闭合处理，常使用三维Smagorinsky方案（Smagorinsky, 1963）或1.5阶三维湍流动能（3-Dimensional Turbulent Kinetic Energy，简称3D-TKE）闭合方案（Deardorff, 1980）来表征次网格湍流的作用。一般而言，模式网格尺度不同，所使用的湍流参数化方案也应不同。次网格湍流混合的参数化方案“本质上”应依赖于模式网格的尺度。

随着目前中尺度数值预报模式的水平分辨率越来越高，当模式网格距Δ与含能湍涡的长度尺度l_e相当时，模式的动力过程可解析一部分湍流运动，而剩余的湍流运动仍需参数化。Wyngaard（2004）将该区域称为“Terra Incognita”（未知领域），通常也被称为湍流参数化的“灰色区域”。次网格混合的参数化方案依赖于模式网格的尺度，“灰色区域”尺度下湍流参数化问题既不同于传统的PBL参数化，也不同于LES参数化。

基于上述问题，人们所期望的是设计一种适应网格距变化的湍流参数化方案。但目前，在“灰色区域”尺度下尚没有较为合适的方案，因此对现有的传统PBL方案进行全面的评估，考察其在“灰色区域”尺度下的表现和适用性及其随水平网格距变化的特征，是改进PBL方案，使其能够适应分辨率变化的前提和基础。在分辨率较粗的中尺度下，国内外学者对PBL方案做了较多评估（邓国等，2005；Hu et al., 2010；邱贵强等，2013；徐慧燕等，2013；张小培和银燕，2013；周强和李国平，2013；黄文彦等，2014；王子谦等，2014; Cohen et al., 2015; 高笃鸣等，2016；Muñoz-Esparza et al., 2016）。但针对PBL方案在“灰色区域”尺度下的评估仍较少：Zhou et al.（2014）设计了一系列对流边界层的理想试验，使用LES及PBL方案（Sun and Chang, 1986）分析了网格距从3600 m向25 m变化过程中的模式表现，结果表明当网格距落入“灰色区域”时，模式动力过程激发了具有尺度依赖性的对流，从而使准确模拟该尺度下的边界层湍流过程变得更加困难。Shin and Dudhia（2016）评估了WRF模式中4种非局地PBL方案[YSU，ACM2（Asymmetric Convective Model 2），EDMF（Eddy Diffusivity Mass Flux），TEMF（Total Energy Mass Flux）]和一种局地PBL方案（MYNN2.5）在次公里网格距下的表现，分析表明以上5种PBL方案均对尺度不敏感。然而，上述工作在PBL方案的选择、所考察的网格距等方面仍不够全面。本研究选取的四种PBL方案，涵盖了PBL参数化常用的两大类方法（K廓线构造方法及湍流动能闭合方法），且四种PBL方案对局地与非局地作用的表征方法各有不同。本研究借助理想干对流边界层的大涡模拟，构造不同分辨率下的参考湍流场，考察了这四种PBL方案在“灰色区域”尺度下的表现和适用性，为发展适用于高分辨率数值模式的次网格湍流混合方案提供了科学参考和实践基础。

2 大涡标准试验与分析方法 2.1 大涡标准试验

为了评估边界层方案在“灰色区域”尺度下的表现，首先需要建立基于高分辨率LES的标准试验，并将其视为模拟大气的“真值”。本研究所选用的LES模式为WRFv3.7的LES版本，使用适当的次网格湍流参数化方案，WRF作为一种LES模式得到了广泛的使用（Moeng et al., 2007; Skamarock et al., 2008）。

本研究利用WRF-LES对理想的干对流边界层进行模拟，模式积分区域为30 km×30 km，水平分辨率为Δ_LES=50 m，顶层高度为2 km，垂直取100层。水平侧边界使用周期性边界条件。试验由固定地表热通量（Q₀=0.24 K m s^-1）和东西向地转风（U_g=10 m s^-1）驱动。LES的次网格湍流参数化采用1.5阶的3D-TKE闭合方案。模式积分4 h，时间步长为0.2 s。初始位温θ廓线为

$ \theta =\left\{ \begin{array}{l} {\rm{300 K, ~~~~~~~~~~~0}} ~~~ z \le 925{\rm{ m, }}\\ 300{\rm{ K}} + \left({z - 925{\rm{ m}}} \right) \times 0.0536{\rm{ K }}{{\rm{m}}^{ - {\rm{1}}}}{\rm{, }}\\ {\rm{ ~~~~~~~~~~~~~~~~~~~~~~~ 925 m}} ~~~ z \le 1{\rm{075 m, }}\\ 308.05{\rm{ K}} + \left({z - 1075{\rm{ m}}} \right) \times 0.003{\rm{ K }}{{\rm{m}}^{ - {\rm{1}}}}, \\ {\rm{ ~~~~~~~~~~~~~~~~~~~~~~~ }}z ~~~ {\rm{1075 m, }} \end{array} \right. $

(1)

其中，z为高度。

2.2 分析方法

基于水平网格距为50 m的LES结果，使用Dorrestijn et al. (2013)的方法构造分辨率从250 m至9000 m的参考湍流场（reference，简称REF），以此作为各分辨率模拟的“真值”。该方法所构造的参考场可提供网格与次网格尺度的湍流统计信息，被广泛用于研究模式在各个分辨率上的湍流特征（Cheng et al., 2010; Honnert et al., 2011; Dorrestijn et al., 2013; Shin and Hong, 2013; Zhou et al., 2014），构造方法如图 1所示：在得到LES的模拟结果后，设LES的水平计算区域为D×D，“灰色区域”内任一水平网格距为Δ，Δ应介于D与Δ_LES之间。此时标准场中包含K个Δ×Δ的次网格，K=(D/Δ)²，而在每个Δ×Δ区域内，有J个Δ_LES×Δ_LES网格，J=(Δ/Δ_LES)²。对于水平网格距为Δ的参考场，每个格点上的变量φ可由对网格内LES结果做平均得到：

$ {\overline \phi ^{{\Delta _k}}}={(J)^{ - 1}}\sum\limits_j {{\phi _{j, k}}} , $

(2)

其中，下标j, k分别代表第k个次网格和该次网格中第j个LES网格。Δ_k表示第k个次网格。

对于任一水平网格距Δ，整个水平计算区域内φ的总通量输送包括网格（可分辨的）与次网格（无法分辨的）部分（下文分别以上标R(Δ)与S(Δ)表示）：

$ \left\langle {w'\phi '} \right\rangle ={\left\langle {w'\phi '} \right\rangle ^{{\rm{R(}}\Delta {\rm{)}}}} + {\left\langle {w'\phi '} \right\rangle ^{{\rm{S(}}\Delta {\rm{)}}}}, $

(3)

其中，w为垂直速度。网格尺度的湍流通量可表示为

$ {\left\langle {w'\phi '} \right\rangle ^{{\rm{R(\Delta)}}}}={K^{ - 1}}\sum\limits_k {\left({{{\overline w }^{{\Delta _k}}} - \left\langle w \right\rangle } \right)\left({{{\overline \phi }^{{\Delta _k}}} - \left\langle \phi \right\rangle } \right)} , $

(4)

次网格尺度通量则由下式表示：

$ \begin{array}{l} {\left\langle {w'\phi '} \right\rangle ^{{\rm{S(\Delta)}}}}=\\ {K^{ - 1}}\sum\limits_k {\left\{ {{J^{ - 1}}\sum\limits_j {\left[ {({w_{j, k}} - {{\overline w }^{{\Delta _k}}})({\phi _{j, k}} - {{\overline \phi }^{{\Delta _k}}}) + f_{j, k}^{w\phi }} \right]} } \right\}} , \end{array} $

(5)

其中，〈〉代表对于整个水平区域D×D的平均。'表示相对于区域平均的扰动。$f_{j, k}^{w\phi }$是Δ_LES的次网格输送量，由LES的次网格参数化方案得到。

用上述方法分别得到Δ=250 m、500 m、1000 m、3000 m、9000 m的参考场，可与PBL方案在各个分辨率下的模拟进行比较。由下文可知，9000 m时次网格热通量占100%，为典型的中尺度分辨率，Δ=3000 m时可解析小部分通量输送，Δ=250~1000 m时参考场的通量输送一部分由次网格参数化表征，一部分可被模式动力过程解析。通过比较这5个网格距上的不同PBL试验结果，可以考察现有PBL方案在“灰色区域”尺度下的表现和适用性。

2.3 边界层参数化方案简介

本文选取了4种有代表性的PBL方案，对其进行评估与对比，分别为YSU方案、MYJ方案、MYNN2.5（Mellor-Yamada-Nakanishi-Niino Level 2.5）和MYNN3方案（Nakanishi and Niino, 2004, 2006, 2009）。具体方案如下：

（1）YSU方案

YSU方案为一阶非局地K闭合方案。相比于局地方案，非局地方案不仅考虑各层之间物质和通量交换，同时也考虑了相邻模式格点上的交换。

对任意变量φ，其通量$\overline {w'\phi '} $可表示为

$ \begin{array}{l} - \overline {w'\phi '} ={K_\phi }\frac{{\partial \Phi }}{{\partial z}} - {K_\phi }b\frac{{{{(\overline {w'\phi '})}_0}}}{{{w_{s0}}{z_i}}} - {(\overline {w'\phi '})_{{z_i}}}{(\frac{z}{{{z_i}}})^3},\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{①}}\;\;\;\;\;\;\;\;\;{\rm{②}}\;\;\;\;\;\;\;{\rm{③}} \end{array} $

(6)

其中，Φ 代表网格尺度变量，φ'代表扰动场，K_φ为湍流扩散系数，z为高度，z_i为边界层高度，${(\overline {w'\phi '})_0}$为地面通量，${(\overline {w'\phi '})_{{z_i}}}$为逆温层通量，w_s0为z=0.5z_i处的混合层特征速度，b为比例系数。等式右侧①项为变量φ的局地梯度项，对于一阶方案，K_φ由经验廓线确定。如其动量扩散系数K_M可表示为

$ {K_{\rm{M}}}=\kappa {w_s}z{\left({1 - \frac{z}{{{z_i}}}} \right)^2}, $

(7)

其中，κ为冯卡尔曼（von Kármán）常数，w_s为混合层特征速度。

公式（6）中①项代表了垂直层次间的湍流交换，②项为非局地的负梯度输送，③项代表了夹卷层湍流输送的作用。②、③两项均为非局地项。为考察局地一阶方案在“灰色区域”下的表现，设计了一组敏感性试验，该组试验对YSU方案进行了调整，使公式（6）中非局地项②项和③项均为0，只保留梯度输送项①，因此命名该组试验为LYSU（Local YSU）。

（2）MYJ方案

MYJ为1.5阶局地TKE（Turbulent Kinetic Energy）方案。与MYNN2.5及MYNN3一样，采用Mellor-Yamada闭合模型（Mellor and Yamada, 1982）。高阶闭合模型因增加了更多的控制方程，从而对湍流输送有更强的物理约束性。

MYJ方案中动量扩散系数K_M的表示方法与K闭合方案不同，即

$ {K_{\rm{M}}}=lq{S_{\rm{M}}}, $

(8)

其中，S_M为稳定度函数，l为湍流混合长，由诊断方程得到：

$ l={l_0}kz/(kz + {l_0}), $

(9)

其中，${l_0}=a\int\limits_0^{{z_i}} {\left({\left| z \right|q{\rm{d}}z} \right)} /\int\limits_0^{{z_i}} {(q{\rm{d}}z)} $，a为常数。

${q^2}/2=\left({1/2} \right)(\overline {u{'^2}} + \overline {v{'^2}} + \overline {w{'^2}})$（u'为东西方向水平风速扰动分量，v'为南北方向水平风速扰动分量，w'为垂直速度扰动分量），为湍流动能，由TKE预报方程得到：

$ \begin{array}{l} \frac{{\partial {q^2}}}{{\partial t}}= - \frac{\partial }{{\partial z}}\overline {w'(u{'^2} + v{'^2} + w{'^2} + 2p'/{\rho _0})} - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{①}}\\ 2(\overline {u'w'} \frac{{\partial U}}{{\partial z}} + \overline {v'w'} \frac{{\partial V}}{{\partial z}}) + 2\frac{g}{{{\Theta _0}}}\overline {w'{\theta _v}'} - 2\varepsilon ,\\ \;\;\;\;\;\;\;\;\;{\rm{②}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{③}}\;\;\;\;{\rm{④}} \end{array} $

(10)

其中ρ₀为空气密度参考态，U为东西方向水平风速网格尺度量，V为南北方向水平风速网格尺度量，g为重力加速度，Θ₀为位温参考态，θ_v为虚位温，ε为耗散率。①项代表TKE的垂直扩散，②项代表TKE的切变产生项，③项代表浮力产生项，④项代表耗散项。因此，这类方案亦被称作TKE方案。本研究中，MYJ、MYNN2.5及MYNN3均为TKE方案。

（3）MYNN2.5方案

MYNN2.5方案同样为1.5阶局地TKE方案。与MYJ方案不同之处主要在于混合长的诊断方程。MYNN2.5方案中混合长的诊断方程较为复杂，如下式：

$ \frac{1}{l}=\frac{1}{{{l_{\rm{S}}}}} + \frac{1}{{{l_{\rm{T}}}}} + \frac{1}{{{l_{\rm{B}}}}}, $

(11)

其中，l由三种不同的长度尺度构成，l_S为近地层长度尺度，l_T为与边界层高度有关的长度尺度，l_B与浮力长度尺度有关，诊断方程表明该混合长于由三种长度尺度中最短的一项所决定。

（4）MYNN3方案

MYNN3方案相比MYNN2.5方案，阶数更高，增加了$\overline {\theta {'^2}} $等高阶量的预报方程，如

$ \frac{{\partial \overline {\theta {'^2}} }}{{\partial t}}= - \frac{\partial }{{\partial z}}\overline {w'\theta {'^2}} - 2\overline {w'\theta '} \frac{{\partial \Theta }}{{\partial z}} - 2{\varepsilon _\theta }, $

(12)

其中，θ为位温，$\overline {w'\theta '} $代表热量通量，ε_θ为$\overline {\theta {'^2}} /2$的耗散率。其中热量通量可表示为

$ \begin{array}{l} - \overline {w'\theta '} =lq{S_{{\rm{H2}}{\rm{.5}}}}\frac{{\partial \Theta }}{{\partial z}} + lq{\Gamma _\theta },\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{①}}\;\;\;\;\;\;\;\;\;\;\;\;{\rm{②}} \end{array} $

(13)

其中，①项代表局地梯度输送，②项代表了负梯度方向的输送，即非局地作用，其中${\mathit{\Gamma } _\theta } \equiv S{'_{\rm{H}}}\partial \Theta /\partial z$，$S{'_{\rm{H}}}$为无量纲稳定度函数。因此亦可认为MYNN3为非局地方案。

综上可知，本研究所选的四种方案涵盖了业务数值模式中常用且计算效率较高的PBL方案，包含了低阶与高阶、局地与非局地方案。具体试验设计如表 1，通过分析比较这四种方案，考察不同类型的PBL方案在“灰色区域”下的表现，并从方案原理出发深入分析与比较“灰色区域”下高阶闭合、非局地项等的作用，为下一步发展“灰色区域”下次网格混合参数化方案提供了参考。

本研究中LES试验的大涡翻滚的特征时间${\tau _*}(\sim {z_i}/{w_*})$约为10分钟，积分4小时后，$t \ge 4{t_0}$（${t_0}=6{\tau _*}$），LES已达到统计意义上的准平衡态（Moeng and Sullivan, 1994）。以Δ=500 m为例，计算该时刻各PBL方案近地层及浮力特征值，如表 2所示。由于不同的边界层方案对应不同的近地层方案，地表摩擦速度u_*不同（表 2），从而导致地表动量通量不同。另外，在干对流边界层下，边界层内湍流运动主要由地表热力强迫引起。因此本研究着重考察热量通量在“灰色区域”下的表现。由于篇幅所限，不在本文中对动量通量做具体讨论。

3 试验结果

本研究中LES试验积分4小时后，已达到统计意义上的准平衡态。因此，选取模式积分4小时后的结果进行比较。

3.1 平均场

图 2为各PBL方案在不同分辨率下平均位温的垂直廓线。其中YSU方案（图 2a）在混合层内的结果与LES（黑色实线）最为吻合，垂直廓线几乎不随水平分辨率变化。其余方案（包括不含非局地项的LYSU）的位温垂直廓线随着网格距的减小逐渐接近LES廓线；夹卷层中，位温随网格距减小而减小；边界层高度随网格距减小而增加。

进一步比较平均位温梯度廓线（图 3），LES的位温梯度在混合层内的垂直分布呈下半层为负，上半层为正，随高度逐渐增加的特点。YSU方案的位温梯度自400 m左右高度向上均为正梯度，即稳定层结；其余方案分辨率较粗时（Δ＞3000 m），混合层内位温梯度均为负，即弱不稳定层结，进入“灰色区域”后（Δ＜3000 m），分布与LES相近。

位温廓线的差异是由各PBL方案对热通量输送的不同表征所造成的。局地方案由于缺少对较大尺度涡动引起的非局地湍流输送的表征，因此必须通过维持负的位温梯度，以保证一定的次网格通量输送。LYSU试验在次公里（Δ≤1000 m）网格距下的位温与位温梯度廓线与TKE方案的形状较为接近，这是由于在“灰色区域”，部分湍流输送，尤其是尺度较大的非局地输送，可被模式动力过程解析，弥补了局地梯度输送的不足。基于TKE的PBL方案（MYJ、MYNN2.5、MYNN3）随分辨率的提高，位温廓线逐渐接近LES。而YSU方案通过非局地项[公式（6）的②③项]表征了较大涡动对通量的输送，使得热量充分混合，混合层内的位温梯度在“灰色区域”内大于0，呈弱稳定层结。

3.2 总通量

当模式水平网格距小于边界层内最大含能湍涡尺度时（Δ＜l_e），总热通量$\left\langle {w'\theta '} \right\rangle _{{\rm{PBL}}}^{{\rm{total}}}$为PBL方案计算的次网格热通量与模式动力过程解析的网格尺度热通量之和：

$ \left\langle {w'\theta '} \right\rangle _{{\rm{PBL}}}^{{\rm{total}}}={\left\langle {w'\theta '} \right\rangle ^{{\rm{R(}}\Delta {\rm{)}}}} + {\left\langle {w'\theta '} \right\rangle ^{{\rm{S(\Delta)}}}}. $

(14)

图 4为各PBL方案在不同分辨率下总热通量的垂直廓线。与位温廓线相似，PBL方案中YSU的廓线几乎不随水平分辨率变化。其余方案的热通量在夹卷层内随水平分辨率有明显的变化，但在混合层内非常接近。

定义边界层高度z_i为总热通量达到最小值的高度，以z_i对纵坐标做归一化处理，计算0.15z_i~0.85 z_i之间PBL方案总热通量$\left\langle {w'\theta '} \right\rangle _{{\rm{PBL}}}^{{\rm{total}}}$相对LES$\left\langle {w'\theta '} \right\rangle _{{\rm{PBL}}}^{{\rm{total}}}$的误差R：

$ R=\frac{{\left\langle {w'\theta '} \right\rangle _{{\rm{PBL}}}^{{\rm{total}}} - \left\langle {w'\theta '} \right\rangle _{{\rm{LES}}}^{{\rm{total}}}}}{{\left\langle {w'\theta '} \right\rangle _{{\rm{LES}}}^{{\rm{total}}}}} \times 100\% . $

(15)

如图 5所示，所有PBL方案的相对误差均较小（小于3%）。其中MYNN2.5方案最接近LES结果，YSU方案次之。Δ=9000 m时，所有方案的误差最大，进入“灰色区域”后，误差明显减小。

上文分析了各PBL方案在“灰色区域”内总通量的变化，下面分别计算0.15z_i~0.85z_i之间（混合层内）网格（${\left\langle {w'\theta '} \right\rangle ^{{\rm{R(\Delta)}}}}$）和次网格（${\left\langle {w'\theta '} \right\rangle ^{{\rm{S(\Delta)}}}}$）热通量与总通量的比值，如图 6所示，点虚线为各PBL方案的结果，实线为REF结果。可以看到，对于网格与次网格通量的分配，各PBL方案差异明显，且对水平网格距敏感。具体表现为，所有PBL方案模拟的次网格通量比例（蓝色点线）均随水平网格距增大而增加，网格通量比例随水平网格距增大而减少。MYNN3方案的网格与次网格通量比例随分辨率的变化与REF最为接近。YSU方案的${\left\langle {w'\theta '} \right\rangle ^{{\rm{S(\Delta)}}}}$比重在“灰色区域”内偏大，当分辨率Δ＞1000 m，${\left\langle {w'\theta '} \right\rangle ^{{\rm{S(\Delta)}}}}$的比例接近100%，几乎全部的${\left\langle {w'\theta '} \right\rangle ^{{\rm{R(\Delta)}}}}$都被抑制。剩余方案则表现为次网格通量比例较REF偏低，网格通量比例偏高。

从图 6 REF的分布可知，Δ=500 m附近，次网格与网格热通量之比约为1:1，可认为该水平网格距为典型的“灰色区域”尺度。以该尺度为例，观察近地层（0.1z_i）位温扰动及垂直速度的水平截面，如图 7所示。从REF的水平分布来看，位温扰动及垂直速度分布介于带状（horizontal rolls）与细胞状（cells）之间。边界层中对流结构取决于浮力作用与机械作用的相对强度ζ，当机械作用占主导时（−25＜ζ＜0）以带状结构为主，当浮力作用主导时（ζ远小于-25）以细胞状结构为主（Deardorff, 1972; Lemone, 1973）。LES中ζ为-28.55（表 2），接近-25，因此带状结构不明显。

对比PBL方案可知，五组试验的ζ均在-25~0之间（表 2），对流结构以带状为主，其中YSU方案与MYNN3方案的ζ较大，分别为-18.27与-18.57，带状对流结构更为明显。YSU方案模拟的对流最弱，MYNN3方案次之，其余三组试验较强，这与它们模拟的网格湍流强弱是相一致的（图 6）。

由此可见，虽然传统的PBL方案在“灰色区域”内能够较为准确地模拟出整个区域的平均位温及总通量，但在细致的结构上，无论垂直分布还是水平分布，仍存在一些问题。在“灰色区域”内，局地与非局地方案，高阶与低阶方案的不同，都会导致次网格与网格通量的不同表现。

3.3 次网格通量

图 8为各PBL方案的次网格通量在不同分辨率下的垂直廓线。可以看到PBL方案模拟的次网格热通量均随分辨率变化，而变化特征各不相同。对比REF廓线可知，Δ=500~3000 m时PBL方案中仅YSU方案模拟出了夹卷层内较明显的负次网格通量，其余方案对该负通量都有显著低估。

在混合层内，PBL方案的次网格热通量均随水平网格距减小而降低，而变化速率各不相同。YSU方案随分辨率的变化幅度最小，在Δ=250~1000 m间YSU方案高估了次网格通量，而当分辨率较粗时（Δ≥3000 m）与REF较为一致；MYNN3方案在Δ≥3000 m下模拟的次网格通量与REF十分接近，至Δ=1000 m时次网格通量明显下降；MYNN2.5方案的次网格通量在Δ=3000~9000 m之间有较明显的变化，该变化特征也可从3.2节的次网格与网格通量比例图上看出（图 6）。与REF相比，每个PBL方案都有各自表现最好的分辨率，如YSU方案在Δ=9000 m和MYNN3方案在Δ=500 m时。

为了进一步揭示PBL方案对分辨率变化的敏感性及其原因，分别考察非局地方案（YSU，MYNN3）中局地项（图 9a、c）和非局地项（图 9b、d）对水平分辨率的敏感性。YSU方案的非局地项包括负梯度和夹卷层输送[公式（6）的②项和③项]，而这些量几乎不随水平网格距变化（图 9b），因此YSU方案的非局地通量对水平网格距不敏感。虽然局地项对水平网格距较为敏感（图 9a），但其所占比例较小，因此总的次网格通量并未表现出较强的分辨率依赖性（图 8a）。对比LYSU试验，去除非局地项后次网格通量表现出对水平网格距较强的依赖性（图 8b）。

MYNN3方案较MYNN2.5方案增加了$\left\langle {\theta {'^2}} \right\rangle $等高阶量的预报方程，在计算热通量时，增加了负梯度输送项[公式（13）②项]，因此混合较MYNN2.5方案更加充分。与YSU方案不同的是，负梯度输送项对水平网格距较为敏感（图 9d），梯度部分输送的次网格通量[公式（13）①项]同样随分辨率的增加而减小（图 9c）。因此，MYNN3方案在“灰色区域”内并没有过强的次网格输送，且对次网格通量的表征具有较好的分辨率依赖性。

从上述分析可知，YSU方案在“灰色区域”内过强的非局地次网格通量导致其对总次网格通量的高估。Shin and Hong（2013, 2015）通过对LES试验结果的诊断分析，定量描述了非局地项与局地项随分辨率变化的特征，并指出次网格通量对分辨率的依赖性主要由非局地项决定，因此在高分辨模式中对非局地项的作用及其随分辨率变化特征的准确描述至关重要。

3.4 网格尺度通量

当模式分辨率逐渐提高，水平网格距逐渐减小，进入“灰色区域”时，部分含能湍涡的尺度超过水平网格距（Δ＜l_e），模式动力过程可部分解析湍流运动，无需再对其参数化。从LES以及REF的可分辨湍流动能(${\left\langle {{\mathop{\rm TKE}\nolimits} } \right\rangle ^{{\rm{R(}}\Delta {\rm{)}}}}$)随时间演变的曲线上（图 10e）可以看出，随着水平网格距增大，可分辨湍流动能减小，Δ≥3000m后，网格距尺度大于最大的含能涡尺度l_e（通常以大气边界层的高度z_i代表最大含能涡的尺度），可分辨湍流动能接近于0。而在Δ=50~1000 m，经过t=t₀（约1小时后），可分辨湍流动能逐渐稳定，之后随地面不断加热而缓慢增加。Shin and Dudhia（2016）认为这表明经过t=t₀后湍流运动进入了准平衡状态。

考察PBL方案的模拟结果发现，四种方案达到准平衡状态的时间均不相同，水平网格距越大，可分辨湍流动能的激发时间越晚，到达平衡的时间也越晚。从图 10a可知，YSU方案模拟出的可分辨湍流动能在所有时次均偏低，Δ≥3000 m后，可分辨湍流动能几乎为0。MYJ方案次之，Δ=3000 m时已接近0。而MYNN2.5方案和MYNN3方案的可分辨湍流动能较为显著。

以上现象可用Rayleigh-Benard对流机制解释。Ching et al.（2014）指出，当水平网格距变小时，模式自身动力过程更易激发“次级环流”（Convectively Induced Secondary Circulations by Model，以下简称M-CISC）。次级环流源自近地面超绝热层热量的积累，当边界层混合不够充分时，近地面累积了过多的热量，当雷利（Rayleigh）数$Ra=g\alpha \beta {d^4}/kv$（其中，α为热膨胀系数，β为垂直温度梯度，d为流体厚度，k为热传导系数，ν为动力粘性系数）达到临界值Ra_c时，M-CISC产生，Ra_c的大小与水平分辨率，即可分辨的最小波长有关。M-CISC的出现提高了热量的输送和混合效率，中和了混合层内原本弱不稳定的静力层结，使得Ra数降低至临界值以下。Zhou et al.（2014）也利用Rayleigh-Benard热力稳定度解释了边界层内动力激发对流的过程。

从上述机制出发，考察M-CISC激发前后边界层内结构的差异，分别对比四种方案在积分1小时后和积分8小时后两个时次的平均位温垂直廓线与垂直速度w的水平分布特征。从图 11的位温廓线结构来看，相比积分1小时后，经过7小时积分后整层位温大幅提升，混合层与夹卷层均增厚。积分1小时后，MYNN2.5方案与MYNN3方案的Δ=3000 m与Δ=9000 m廓线在混合层上半层均表现为弱不稳定层结。至积分8小时后，由于M-CISC发展旺盛，Δ=3000 m在混合层上半层转为稳定层结，Δ=9000 m因缺乏M-CISC对湍流输送的补充，位温廓线仍呈不稳定层结。

从近地面层（0.1z_i）垂直速度的水平分布来看（图 12），至积分8小时后，Δ=250~3000 m试验的对流强度明显强于积分1小时后，尤其是MYNN2.5方案与MYNN3方案的Δ=3000 m试验，网格上的对流发展旺盛。此外，对流结构的尺度变大。而Δ=9000 m试验的对流强度没有变化，因为在该分辨率下的Ra_c在积分8小时后仍未达到。

另外，随着时间增加，若满足了M-CISC的触发条件，可分辨热通量增加，而次网格热通量所占比例受到相应抑制。从可分辨热通量与次网格热通量分配（图略）来看，Δ=3000 m下时各PBL方案${\left\langle {w'\theta '} \right\rangle ^{{\rm{R(\Delta)}}}}$均在7小时后明显增大。相应地，“灰色区域”下次网格通量比例减小，表明其受到了过强的可分辨湍流输送的抑制。从3.2节分析可知，网格尺度通量输送与次网格通量互为补充，例如，当PBL方案高估了次网格通量时，网格尺度通量得到低估，从而保持总通量的一致（图 4），反之亦然。因此，由于不合理的M-CISC产生了过强的可分辨热通量，抑制了次网格通量，从而影响PBL方案对次网格通量的表征。

相较TKE方案，YSU方案中激发的M-CISC强度较弱（图 11），在较小的网格距尺度下（Δ ≤ 1000 m）才开始出现较显著的可分辨湍流动能（图 10）。这可能是由于该方案包含了非局地输送，次网格混合足够充分，近地面热通量难以积累，无法达到Ra_c，从而无法触发M-CISC，较其他方案不易激发旺盛的对流。为考察不同强度的次网格混合对可分辨湍流动能的抑制效果，利用YSU方案设计敏感性试验如表 3所示。

如表 3所示，设计了3组敏感性试验，将YSU方案中非局地项中的夹卷项（公式（6）③项）关闭，令输送近地面热量及动量的负梯度项（公式（6）②项）（Counter-gradient Term, 简称CT）为控制试验的0倍、1倍与2倍（0CT、1CT、2CT试验）。结果如图 13所示，控制试验CTRL与1CT的差别只来自于夹卷层造成的非局地输送，因此可分辨湍流动能曲线较为接近。非局地输送为0（0CT）时，可分辨湍流动能在第一个小时增至约0.4 m² s^-2，之后缓慢上升。当负梯度项倍数增加，非局地混合增强，可分辨湍流动能减小。在2CT试验中，混合过于强盛，可分辨湍流动能降至0，在8小时内被完全抑制。可见，过强的非局地次网格通量，不但抑制了局地次网格通量（3.3节），而且抑制了网格通量的产生。

由上文可知，进入“灰色区域”，随着水平网格距的减小，Ra_c降低，更易激发M-CISC。水平网格距决定了Ra_c的大小，初始时刻次网格湍流输送的强度决定了近地面层滞留热量的多少，由此可确定Rayleigh-Benard对流激发的条件和时间。可见，可分辨湍流输送的强度与时间演变一方面与水平网格距有关，另一方面与不同的参数化方案选取有关。增大网格距与近地面热量的输送效率，都可抑制M-CISC的产生。然而次级环流与中小尺度天气密切相关，随着业务预报模式分辨率越来越高，需要准确表达“灰色区域”内可分辨湍流的特征。

对比四种PBL方案在次网格湍流与网格湍流的表现可知，不同的参数化方法，首先直接影响了次网格通量的模拟效果。所有方案模拟的混合层内次网格热通量相对总通量的比例都随分辨率提高而减小，与REF的变化趋势一致，只是减小的幅度和速率各有不同。同时，PBL方案还决定了边界层发展初期的湍流输送与混合强度。由图 6可知，局地方案在“灰色区域”内低估了次网格湍流输送，不充足的输送导致近地面滞留了过多的热量，从而提供了动力触发对流的条件。这种由动力过程触发的对流，不仅会掩盖真实对流，如因下垫面不同导致的次级环流（Ching et al., 2014），还会抑制次网格湍流输送。并且动力触发对流的条件随分辨率提高而降低，其触发时间随分辨率提高而提前，增加了数值模式准确表征边界层湍流过程的难度。

3.5 PBL方案适用性评估

由前文可知，PBL方案在“灰色区域”尺度下的适用性是与具体尺度有关的。由图 2可知，在“灰色区域”内，YSU的位温廓线始终与LES吻合较好，其余方案随着网格距的减小也逐渐接近LES廓线；图 6与图 8则表明，混合层内，对于每种方案，均存在一个网格距，在该网格距下次网格及网格热通量分配与LES最为接近。因Δ=500 m时次网格与网格热通量之比约为1:1（见图 6），为典型的“灰色区域”尺度，本节首先以Δ=500 m为例，比较该分辨率下各PBL方案的适用性。由图 10可知，模式积分4小时时，Δ=500 m下各组试验受M-CISCs影响不显著，且达到了LES的准平衡态，因此下文选取模式积分4小时的资料进行比较。

图 14为各PBL方案在Δ=500 m时平均位温的垂直廓线。从图中可知，在混合层中，YSU与LES最为吻合，MYNN2.5与LYSU次之，表现为整层偏冷约0.2 K，MYJ与MYNN3模拟的混合层位温廓线较LES相差最大，约为0.4 K。

图 15a为Δ=500 m时各PBL方案总热通量的垂直廓线。在混合层内5种PBL方案与REF吻合均较好，但夹卷层内差异较大。结合表 2可知，五种PBL方案模拟出的边界层高度z_i（总热通量达到最小值的高度）均低于REF。YSU最为接近，之后为MYNN2.5、LYSU、MYJ及MYNN3。从夹卷率R_E来看，YSU与REF吻合较好，其余方案均偏弱，MYJ最弱。

分析次网格热通量廓线（图 15b）可知，在700 m高度以下，MYNN3与REF廓线最为接近，YSU次网格通量偏高，其余三种方案偏低。在夹卷层内，仅YSU模拟出了夹卷层内较明显的次网格负通量，MYNN3次之，其余三种方案均显著低估了该部分次网格负通量。对比网格热通量廓线（图 15c），MYNN3在700 m高度以下与REF廓线最为接近，YSU方案偏弱，其余三种方案偏高。

由上述分析可知，从位温与总通量廓线来看，Δ=500 m时YSU的模拟最为接近LES，但对次网格及网格通量的分配欠佳；在混合层中下层MYNN3对次网格和网格通量的表征较好，混合层中上层至夹卷层各PBL方案均不理想。即这些方案中不存在一种最佳方案，能够既准确表征边界层的热力结构，又给出最接近REF的湍流统计特征。

类似地，当分辨率改变时，如提高至Δ=250 m，此时MYNN2.5对位温结构的模拟更为准确，而在次网格与网格热通量表征上欠佳；MYJ在混合层内湍流统计特征最为接近LES，而其他方面则不够准确。在“灰色区域”尺度下，本文所考察的PBL方案中不存在一个表征热力结构与湍流统计特征均最佳的方案。这是因为传统PBL方案是针对分辨率较粗的中尺度模式而设计的，当初并未考虑分辨率更精细的“灰色区域”尺度下的适用性问题。

近年来一些工作开始尝试改善传统PBL方案在“灰色区域”下的表现。Shin and Hong（2015）基于LES理想实验确定了具有分辨率依赖性的函数，分别对YSU局地和非局地项进行了尺度依赖性的调整；Ito et al.（2015）通过在特征长度表达式中增加了具有尺度依赖性的函数，对Mellor-Yamada 3阶模型做了改进。可见，引入尺度依赖性函数，是现阶段改进传统PBL方案的主要方法。然而目前尺度依赖性函数的建立是基于LES模拟的干对流边界层个例，对于情况更为复杂的真实大气，这些函数是否适用，还有待进一步研究。

4 总结与讨论

本研究利用WRF LES版本对理想的干对流边界层进行高分辨率模拟，基于LES结果构造水平分辨率从250 m至9000 m的参考湍流场，以此作为“真值”，比较了四种业务常用的PBL方案在“灰色区域”下的表现。结果表明：

混合层内总热通量对所使用的参数化方案和水平分辨率均不敏感。不同参数化方案中次网格与网格通量的比例表现出对水平网格距不同的依赖性。网格尺度通量输送与次网格通量互为补充，PBL方案高估了次网格通量时，网格尺度通量被低估，从而保持总通量的一致，反之亦然。

局地PBL方案（LYSU、MYJ、MYNN2.5）在混合层内的平均位温随网格距减小而增大，次网格通量随网格距减小而减小。较参考湍流场，局地方案对次网格通量有所低估，表现为热量向上输送不足，混合不充分。

对两种非局地PBL方案随分辨率变化的敏感性及其原因的分析表明：YSU方案的非局地项几乎不随水平网格距变化而变化，因此该方案对次网格通量的表征并未表现出较强的分辨率依赖性；MYNN3方案的非局地次网格通量（负梯度输送项）随分辨率的增加而减小，对次网格通量的表征具有较好的分辨率依赖性。

YSU方案过强的非局地次网格输送使混合层内温度层结呈弱稳定，抑制了可分辨湍流输送，不易于激发次级环流。

PBL方案在“灰色区域”尺度下的适用性与具体分辨率有关。以Δ=500 m为例，四种PBL方案中不存在一种最佳方案，能对边界层的热力结构和湍流统计特征均有准确的描述。

本研究讨论的PBL方案均针对分辨率较粗的中尺度数值预报模式和大气环流模式而设计，在较粗的分辨率下，可认为大气边界层内所有的湍流运动均为次网格尺度，需对所有的湍流进行参数化，因此不存在湍流参数化“灰色区域”的问题。当模式分辨率逐渐提高时，模式动力过程可分辨部分湍流运动，从而导致了“灰色区域”问题的产生。另外，随着数值天气预报模式的分辨率越来越高，能否准确预报次级环流也是评判模式性能的重要指标之一，此时水平均一性假设不再成立，水平方向的湍流混合过程不应忽略，只考虑了垂直混合的一维PBL方案将不再适用。为解决在“灰色区域”下湍流参数化的问题，可从两个方面入手。一方面可在PBL方案中引入具有分辨率依赖性的函数来改善次网格通量与网格通量的比例；另一方面需在PBL方案中考虑三维混合的作用。

本研究针对干对流边界层，初步分析了一维PBL方案在“灰色区域”尺度下的表现，没有考虑水汽在混合过程中的作用。湿过程中，水汽相变可引发一系列复杂物理过程，参与并影响湍流混合，从而增加边界层模拟的复杂性。因此，高分辨率数值预报模式中PBL方案与其他物理过程的相互作用值得进一步深入研究。