引言

Stull^[1]将大气边界层定义为直接受下垫面影响的那部分大气，它在自由大气和下垫面的动力和热力相互作用中起着十分重要的作用。太阳辐射是大气运动的主要动力来源，下垫面吸收短波辐射后，下垫面和自由大气之间的热量交换主要通过边界层的湍流热量输送完成。边界层内的湍流动量输送能有效地耗散大气运动，从而不会使其无限制地发展。边界层内的物质和水汽输送对云和降水的形成，污染物的扩散以及大气化学变化都有着非常重要的作用。此外，边界层内的低云可有效地反射太阳短波辐射，对全球辐射收支平衡也有着重要的作用。鉴于大气边界层在天气和气候系统中的重要作用，在数值模式中对边界层内物理过程的合理描述是十分必要的，但由于模式本身的垂直和水平分辨率限制和边界层内湍流运动的复杂性，边界层内的一些物理过程不能被模式中的动力框架描述出来，因此需要将这些物理过程进行参数化，即将一些次网格尺度的变量用网格可解析的变量描述出来，使模式的模拟结果中能够体现这些过程的影响。合理的边界层参数化方案在一个成熟的数值模式中是不可或缺的。

数值模式中通过雷诺平均的形式将风速、水汽、位温等物理量分解为网格尺度的平均量和次网格尺度的扰动量，再将其代入原始的纳维-斯托克斯方程组，方程组会因为多了湍流通量项而不能闭合。在边界层中垂直方向湍流输运的重要性要远大于水平方向，边界层参数化方案主要解决的就是垂直湍流通量项的闭合问题，该问题可用公式(1)简单地描述出来：

$ \frac{\partial \overline{\chi}}{\partial t}=-\frac{\partial \overline{w^{\prime} \chi^{\prime}}}{\partial z} $

(1)

其中$\overline{\chi} $表示模式网格点上的物理量，而$\overline{w^{\prime} \chi^{\prime}}$即为需进行参数化的垂直湍流通量项。

从闭合框架的角度来划分，常用的边界层参数化方案可分为1阶闭合，1.5阶闭合以及2阶或更高阶闭合。闭合阶数一般是通过引入的最高阶预报方程来确定的，1阶闭合直接将湍流通量项描述为已知量的函数，1.5阶闭合在描述湍流通量项时引入了部分二阶矩(例如湍流动能)的预报方程，2阶闭合引入了全部二阶矩的预报方程，更高阶的闭合可以此类推。

影响边界层内垂直湍流通量项的物理过程主要有下垫面的湍流输运、边界层顶的夹卷以及边界层顶层积云的长波辐射冷却等。浮力和底层风切变对下垫面的湍流输运有着至关重要的影响，依据浮力对湍流混合是增强还是抑制作用可将大气边界层划分为不稳定边界层(也称对流边界层)和稳定边界层。不稳定边界层多发生在晴朗的白天，下垫面受太阳辐射加热，浮力和风切变共同促进湍流混合；而稳定边界层多发生在夜间，下垫面自然冷却，仅有风切变促进湍流混合，浮力反而对湍流混合具有抑制作用。这两类边界层具有非常不同的内部结构，在一些参数化方案中也进行了区别对待。图 1表示不稳定边界层和稳定边界层的内部结构^[2]。在稳定边界层中，湍流主要由风切变产生和维持，具有较强的局地性，不能形成较大尺度的涡旋；在不稳定边界层中，由于浮力对湍流的促进作用会出现许多上冲的热泡，这些热泡在垂直方向上具有较大的尺度，有些甚至可以跨越整个边界层，它们将近地面的空气带到边界层的中上部，使边界层内部充分混合，因此不稳定边界层内的湍流输送具有很强的非局地性，很多工作对边界层参数化方案中如何考虑这种非局地性进行了研究^[3-6]。在不稳定边界层的顶部，上冲热泡由于惯性作用可突破顶部逆温的限制进入到自由大气，通过周围空气的下曳作用将自由大气中的空气卷入边界层中，使其与边界层内的空气进行混合，这个过程被称为不稳定边界层顶部的夹卷过程。夹卷过程对不稳定边界层的发展有十分重要的作用，对该过程的深入研究一方面是通过对实际观测资料的分析^[7-9]，另一方面通过实验室的理想实验和大涡模拟等数值模拟技术^[10-14]。根据边界层顶是否有云覆盖可将边界层分为晴空边界层和有云覆盖的边界层，特别是有海洋层积云覆盖的边界层。海洋层积云主要存在于副热带大陆西岸的冷洋面之上，顶部具有较强的长波辐射冷却^[15]，这也是它区别于其他种类云的一个显著特征。层积云顶的长波辐射冷却过程会产生自上而下的湍流混合，从而成为除下垫面之外的另外一个驱动边界层内湍流维持和发展的重要来源^[16-17]。Wood^[18]对海洋层积云的物理性质及其与边界层过程之间的相互作用进行了详细的综述。

本文将从闭合框架以及对边界层内重要物理过程的考虑这两个层面总结回顾常用的边界层参数化方案，介绍了几个主要数值预报模式中的业务方案，并且对边界层方案的发展方向进行简要概述。

1 闭合框架 1.1 1阶闭合方案

1阶闭合方案多采用K理论，假设湍流运动与分子运动具有一致性，得出湍流通量是湍流输运系数与网格尺度物理量梯度的乘积，以纬向风速u和位温θ为例：

$ \overline{w^{\prime} u^{\prime}}=-K_{m} \frac{\partial \overline{u}}{\partial z} $

(2)

$ \overline{w^{\prime} \theta^{\prime}}=-K_{h} \frac{\partial \overline{\theta}}{\partial z} $

(3)

其中K_m和K_h分别代表动量和热量的湍流输运系数。

1.1.1 K(Ri)方案

此类方案采用经典的普朗特半经验公式，又被称为混合长方案：

$ K_{m, h}=l^{2}\left[\left(\frac{\partial \overline{u}}{\partial z}\right)^{2}+\left(\frac{\partial \overline{v}}{\partial z}\right)^{2}\right]^{1 / 2} F_{m, h}(R i) $

(4)

l表示湍流混合长，F_{m, h}(Ri)为理查森数Ri的函数，Ri是表征边界层内稳定度的一个重要的无量纲数：

$\ R i=\frac{g}{\overline{\theta}_{v}} \frac{\partial \overline{\theta}_{v} / \partial z}{(\partial \overline{u} / \partial z)^{2}+(\partial \overline{v}/\partial z)^{2}} $

(5)

Louis^[19]和Louis et al.^[20]利用近地层的莫宁-奥布霍夫相似律^[21-24]导出了函数F_m,h(Ri)并将其扩展到了整个边界层, Beljaars^[25]对稳定边界层下的函数F_{m, h}(Ri)进行了改进。类似莫宁-奥布霍夫相似律的形式，Nieuwstadt^[26]提出了局地相似律，Duynkerke^[27]利用该理论导出了F_m,h(Ri)。此外，混合长l多采用Blackadar^[28]的形式：

$ \frac{1}{l}=\frac{1}{\kappa z}+\frac{1}{l_{0}} $

(6)

其中κ为卡曼常数，l₀为渐进混合长，根据需要可取常数值或函数形式。

1.1.2 K廓线方案

此类方案直接预先给定湍流输运系数K的垂直廓线形式^{[5, 29]}：

$ K_{m, h}=\kappa u_{*}\phi_{m, h}^{-1} z\left(1-\frac{z}{h}\right)^{2} $

(7)

其中u_*表示摩擦速度，h表示边界层高度，ϕ_m和ϕ_h是莫宁-奥布霍夫相似律中关于z/L的函数，L表示莫宁-奥布霍夫长度。不同于K(Ri)方案，该类方案在计算湍流输运系数时引入了边界层高度，Seidel et al.^[30]利用观测、模式和再分析资料对多种计算边界层高度的公式进行了检验，并指出Vogelezang and Holtslag^[31]提出的通过整体理查森数诊断得到边界层高度的方法能较好地适用于各类情况。

1.2 1.5阶闭合方案

1.5阶闭合方案在数值模式中应用最广泛的是K(TKE)方案^[32-34], K(TKE, ε)方案虽然也有应用^[35]，但并不多见，因此该小节主要对K(TKE)方案进行介绍。K(TKE)方案仍采用K理论(见公式(2)和(3))，但是在计算湍流输运系数的过程中引入了湍流动能e的预报方程：

$ \begin{align} \frac{\text{d}e}{\text{d}n}=\underbrace{-\overline{u'w'}-\frac{\partial \overline{u}}{\partial z}-\overline{u'w'}\frac{\partial \overline{v}}{\partial z}}_{S}+\underbrace{\frac{g}{{{\overline{\theta }}_{v}}}\overline{w'\theta {{'}_{v}}}}_{B}- \\ \ \ \ \ \ \ \ \ \underbrace{\frac{\partial }{\partial z}\left( \overline{w'e}+\frac{1}{\overline{\rho }}\overline{w'p'} \right)}_{T}-\varepsilon \\ \end{align} $

(8)

湍流动能e是表征湍流强度的重要物理量，S表示湍流的切变产生项，B表示湍流的浮力产生项，T表示湍流输运项和气压相关项，ε表示湍流耗散项。在K(TKE)方案中湍流输运系数K的公式可表示为：

$ K_{m, h}=c l e^{1 / 2} $

(9)

其中c为一个常数，l表示湍流混合长。不同的方案针对混合长l以及湍流动能方程(公式(8))中的各项有不同的假设^[32-34]，但基本框架均大致相同。

1.3 2阶闭合方案

这类方案引入了所有二阶矩的预报方程，并且为了使这些方程闭合，通过一定的假设和简化，将方程中出现的三阶矩转换成由二阶矩和平均量表示的形式，然后联立这些预报方程计算出湍流通量。Mellor and Yamada^[36-37]建立了完整的2阶闭合的参数化框架，并通过进行不同的假设得到不同复杂度的参数化方案。采用他们文章中的术语^[37]，依据不同的假设将方案复杂度分为Level 4、Level 3、Level 2.5和Level 2等四个级别。Level 4的复杂度需求解完整的关于二阶矩的预报方程组，Level 3的复杂度仅需求解湍流动能e和位温扰动的方差$\overline{\theta^{\prime 2}}$的预报方程，Level 2.5的复杂度仅需求解湍流动能e的预报方程，此时方案变为K(TKE)方案，但是计算湍流输运系数的公式变为：

$ K_{m, h}=l e^{1 / 2} S_{m, h} $

(10)

其中S_m，h为$G_{m}=\frac{l^{2}}{e}\left[\left(\frac{\partial \overline{u}}{\partial z}\right)^{2}+\left(\frac{\partial \overline{v}}{\partial z}\right)^{2}\right]$和${{G}_{h}}=\frac{{{l}^{2}}}{e}\frac{g}{{{\overline{\theta }}_{v}}}\frac{\partial {{\overline{\theta }}_{v}}}{\partial z}$的函数，Level 2的复杂度则使方案变为1阶闭合方案中的K(Ri)方案，并且F_{m, h}(Ri)的形式可以直接从对2阶闭合模型所做的假设与简化中推导出来。后续有许多工作对上述文章中的公式和假设进行了改进^[38-40]，Galperin et al.^[38]利用湍流动能的准平衡模型，针对Level 2.5中S_m，h的形式做了简化，使其仅是G_h的函数，Grenier and Bretherton^[41]，Bretherton and Park^[42]在此基础上利用公式(10)对边界层进行参数化；Janjić^[39]通过给混合长加一个上限解决Level 2.5中出现的奇异性问题，并且在Level 2.5的具体计算步骤上兼顾了计算精度与计算效率；Nakanishi and Niino^[40]总结了他们在Level 3的基础上进行的改进，主要有改进了混合长的诊断公式^[43]，考虑了湿过程^[44-45]及提高了Level 3在数值模式中的稳定性^[46]等方面，并且在改进的Level 3基础上也有了改进的Level 2.5方案。

2 重要物理过程 2.1 湿过程

相对于晴空边界层，有云覆盖的边界层因为要考虑边界层与云的相互作用，所以需要使用湿绝热守恒量q_t和θ_l来替代干绝热守恒量q和θ，若模式中使用的是干静力能s=c_pT+gz而不是θ，则用s_l替代。这样做的好处是可以在描述湍流输运的过程中不需要显式地表达出蒸发和凝结过程^[47]，湿绝热守恒量可以表示为：

$ q_{t}=q+q_{l} ; \theta_{l}=\theta\left(1-\frac{L_{v} q_{l}}{c_{p} T}\right) ; s_{l}=c_{p} T+g z-L_{v} q_{l} $

(11)

其中q_l表示液态水混合比，L_v为水汽蒸发潜热系数，c_p为大气比热。若边界层方案中需要湍流通量$\overline{w^{\prime} \theta_{v}^{\prime}}$的计算(例如公式(8)中的B)，则需要将$\overline{w^{\prime} \theta_{v}^{\prime}}$改写成由$\overline{w^{\prime} \theta_{l}^{\prime}}$和$\overline{w^{\prime} q^{\prime}_{t}}$表示的形式^{[41, 44-45, 48-50]}；若方案中使用的是虚干静力能s_v=c_pT_v+gz，而不是虚位温θ_v，需要将$\overline{w^{\prime} s^{\prime}_{v}}$改写成由$\overline{w^{\prime} s^{\prime}_{l}}$和$\overline{w^{\prime} q^{\prime}}$表示的形式^[15，42]。此外，若方案中用到理查森数(公式(5))，则需将其改写为由q_t和θ_l(或s_l)表示的湿理查森数^{[34, 42, 51]}。

2.2 不稳定边界层的非局地性

本文引言中提到不稳定边界层内的湍流输运具有很强的非局地性，但是针对参数化方案中经常采用的K理论(公式(2)和(3))，湍流输运仅与物理量的局地垂直梯度有关，并且都是沿着梯度的方向输运，这是不符合实际的。由图 2可以看出由于大尺度湍涡导致的非局地输运，在不稳定边界层的中上部可能存在逆梯度输运的情况。为了考虑这种非局地性质，K理论的公式可改写为^{[4-5, 32-33]}：

$ \overline{{{w}^{\prime }}{{\chi }^{\prime }}}=-{{K}_{\chi }}\left( \frac{\partial \bar{\chi }}{\partial z}-{{\gamma }_{\chi }} \right) $

(12)

其中γ_χ为逆梯度项，与对应物理量的下垫面湍流通量有关。

除了在K理论的公式中加入逆梯度项这一方法之外，方案中还可采用直接混合的方法来描述不稳定边界层的非局地特性，Pleim^[6]的方案就是这种方法的一个典型例子。如图 3所示，方案假设对流输送源于模式的最低层，从最低层向上直接输送到其他每一层，并且用逐层的向下输送表示补偿下沉运动，方案还用K理论来表示局地的湍流输送：

$ \begin{align} \frac{\partial {{{\bar{\chi }}}_{i}}}{\partial t}=\underbrace{Mu{{{\bar{\chi }}}_{1}}}_{\text{非局地输送}}+\underbrace{M{{d}_{i+1}}{{{\bar{\chi }}}_{i+1}}\frac{\Delta {{z}_{i+1}}}{\Delta {{z}_{i}}}-M{{d}_{i}}{{{\bar{\chi }}}_{i}}}_{\text{补偿下沉}} +\\ \underbrace{\frac{1}{\Delta {{z}_{i}}}\left[ \frac{{{K}_{i+1/2}}\left( {{{\bar{\chi }}}_{i+1}}-{{{\bar{\chi }}}_{i}} \right)}{\Delta {{z}_{i+1/2}}}-\frac{{{K}_{i-1/2}}\left( {{{\bar{\chi }}}_{i}}-{{{\bar{\chi }}}_{i-1}} \right)}{\Delta {{z}_{i-1/2}}} \right]}_{\text{局地输送}} \\ \end{align} $

(13)

其中下标代表模式层，Mu和Md表示混合系数。

另外一种表示不稳定边界层非局地特性的方法借鉴了浅对流方案中的质量通量方案^[52]，在不稳定边界层的情形下将局地的K理论与非局地的质量通量结合起来^[53-58]：

$ \overline{w^{\prime} \chi^{\prime}}=-K_{\chi} \frac{\partial \overline{\chi}}{\partial z}+M\left(\chi_{u}-\overline{\chi}\right) $

(14)

其中M表示质量通量，下标u表示边界层内的强上升气流，如图 4所示。这类方案统称为EDMF(Eddy-Diffusivity Mass-Flux)方案。

2.3 层积云顶的长波辐射冷却和边界层顶的夹卷过程

本文引言中提到层积云顶的长波辐射冷却和不稳定边界层顶的夹卷过程对于边界层的发展有十分重要的作用。对于有层积云覆盖的边界层，云顶辐射冷却是驱动边界层内湍流发展的一个重要来源，在边界层方案中主要有两种方式表示这个过程。在1阶闭合的框架下一般采取双K廓线的方法，即除了预先给定类似公式(7)所示的这条K廓线外，在其基础上再叠加一条K廓线，这条K廓线表示云顶辐射冷却作为湍流动力来源在边界层内产生的湍流混合^[59-62]，其中云顶辐射冷却的强度在模式中通常由辐射通量的散度表示。在1.5阶、2阶和高阶闭合的框架下，方案中通常将云顶辐射冷却产生的浮力项直接写入湍流动能和高阶矩的收支方程中，在方程中将该过程隐式地表达出来^{[42, 63]}。

针对不稳定边界层顶的夹卷过程，通用的参数化框架可表示为如下形式：

$ \overline{\left(w^{\prime} \chi^{\prime}\right)}_{h}=-w_{e} \Delta \overline{\chi} $

(15)

其中$\overline{\left(w^{\prime} \chi^{\prime}\right)}_{h}$代表边界层顶的夹卷通量，Δχ代表物理量在边界层顶之上和边界层顶之下的两个模式层上的差值，w_e代表夹卷速度，针对夹卷过程的参数化可转化为对夹卷速度的计算。多数边界层方案在计算湍流通量的过程中将夹卷通量转化成了对应的湍流输运系数，由公式(2)、(3)和公式(15)可以看出夹卷通量对应的湍流输运系数K_e可以表示为K_e=w_eΔz，也有一些方案在计算湍流通量的过程中直接显式地引入了夹卷通量^[62，64-65]。许多研究^{[42，64，66-73]}依据观测或大涡模拟的结果对夹卷速度提出了相应的计算公式。尽管这些公式的形式差别很大，但是总体来说，夹卷速度正比于湍流能量项(例如下垫面湍流热量和动量通量、边界层内湍流浮力项的积分、边界层内平均湍流动能、云顶长波辐射冷却、云水蒸发冷却等)，反比于边界层的特征长度(例如边界层高度、边界层内的湍流混合长等)以及边界层顶附近的稳定度(例如边界层顶上下的虚位温或虚静力能的差值等)。

3 主要数值模式中的业务方案

在早期的数值模式中大多参照公式(2)和公式(3)所示的K理论，使用计算较为简单和稳定的1阶闭合方案，后续又有很多模式通过使用1.5阶的闭合方案，引入湍流动能这一重要物理量，更细致地对边界层进行描述。但是传统的K理论方案属于局地方案，并不能很好地模拟出不稳定边界层的非局地性。为避免此问题，二阶或更高阶的闭合方案得以在模式中开始进行应用，这类方案通过增加多个湍流预报方程能较真实地反映边界层的特征，但是这类方案计算量较大，对多数模式特别是全球模式来说，计算过于复杂，同时引入了很多的不确定参数，对模式的调试也带来了很大的挑战。在此之后，越来越多的非局地方案在模式中得到应用和发展，它在克服传统的局地方案缺点的同时，也不增加太多的计算量(见本文2.2部分)。

数值模式是一个复杂的系统性工程，在一个模式中使用哪种边界层方案要根据模式的动力框架、模式中其他各类参数化方案、模式分辨率、模式具体模拟结果等多种因素综合考虑，并没有哪一种方案可以适应于所有模式。但是综合来说，使用1阶或1.5阶闭合的非局地方案，并在方案中引入不稳定边界层顶的夹卷和层积云顶的长波辐射冷却，是目前较为主流的形式。另外随着计算能力的迅速发展，2阶或更高阶的方案也在模式中得到更为广泛的应用。

本文选取了在日常天气预报业务中预报员使用较多并且整体预报性能较好的三个业务模式，介绍它们所使用的边界层方案。

3.1 GRAPES_Meso

GRAPES是中国气象局自主开发建立的多尺度通用资料同化与数值预报系统，其中GRAPES_Meso模式是GRAPES的区域中尺度预报模式，该模式的水平分辨率为0.1°，垂直方向有50层^[74]。它使用的边界层方案主要有三个^[75]，分别是如公式(7)和公式(12)所示的使用K廓线并在不稳定边界层中加入了逆梯度项的MRF方案^[5]，在MRF方案的基础上考虑不稳定边界层顶部夹卷的YSU方案^[64-65]以及如公式(10)所示的使用TKE进行1.5阶闭合的MYJ方案^[39]。

3.2 ECMWF IFS-HRES

ECMWF IFS是欧洲中期天气预报中心研发的一套集成预报系统，HRES是系统中的高分辨率大气模式，该模式的水平分辨率约为9 km，垂直方向有137层。如图 5所示，模式根据边界层的性质不同使用不同的方案^[76]。混合层表示不稳定边界层中充分混合的部分，混合层内使用的边界层方案属于如公式(14)所示的EDMF方案。若边界层内无云或者顶部被层积云覆盖，混合层高度即为边界层高度，混合层中的大尺度非局地混合通过EDMF方案中的质量通量项描述；若边界层内存在浅积云，混合层高度为浅积云的云底高度，此时EDMF方案中的质量通量项变为零，非局地混合由浅积云中的浅对流方案描述。混合层中的局地混合通过湍流输运系数K来描述，K的计算方案属于如公式(7)所示的K廓线方案，并且方案还考虑了混合层顶部的夹卷过程以及通过双K廓线的形式考虑了层积云顶的长波辐射冷却过程。稳定边界层以及不稳定边界层中混合层之上的部分使用的方案均属于如公式(4)所示的K(Ri)方案。

3.3 NCEP GFS

NCEP GFS是美国国家环境预报中心研发的全球预报系统，其中它的高分辨率版本的大气模式水平分辨率约为13 km，垂直方向有64层。模式根据边界层的稳定度不同使用不同的边界层方案^[77]，模式将边界层分为强不稳定边界层、弱不稳定边界层、弱稳定边界层以及强稳定边界层。强不稳定边界层使用的方案属于如公式(14)所示的EDMF方案，方案的局地混合部分由公式(7)所示的K廓线方案描述；弱不稳定边界层使用的方案属于K廓线方案，同时加入了如公式(12)所示的逆梯度项；弱稳定边界层使用的方案也属于K廓线方案；强稳定边界层使用的方案属于如公式(4)所示的K(Ri)方案。

4 发展方向

尽管边界层方案历经多年发展，各方面已考虑得较为完善，但是随着模式水平的不断发展，仍有需要改进和继续深入研究的方向。首先是对边界层(特别是有云覆盖的边界层)内各物理过程的更为深入和细致地研究，深化对内在物理机制的认识。其次是考虑边界层方案和其他物理参数化方案(特别是对流参数化方案和云物理方案)之间的联系，建立一个统一的参数化框架，使得模式在边界层湍流输送，浅对流和深对流发展，以及云的形成等各主要物理过程之间能够更加自然地切换。前文所述的EDMF方案就为建立统一的参数化框架提供了很好的基础^[57-58]，另外也有研究在高阶闭合的框架下将各物理方案统一起来，CLUBB(Cloud Layers Unified by Binormals)方案^[78-81]在3阶闭合的框架下以多元高斯分布为基础进行闭合假设，通过求解八个二阶矩和一个三阶矩的预报方程组将边界层、浅对流以及云物理方案统一到同一个计算框架内。此外，随着模式分辨率的不断提高，当模式的网格距小于1 km时，模式的动力框架便可部分分辨湍流过程，这时传统的边界层方案可能就不太适用了，如何处理这种高分辨率下边界层的参数化也是未来研究的主要方向之一^[82-84]。