1 引言

边界层高度作为表征大气边界层的一个重要参数, 广泛用于大气模式边界层参数化、边界层结构以及物理量廓线的研究工作中 (惠小英等, 2011); 大气边界层高度决定了污染物垂直混合的范围, 直接关系到区域内的污染物扩散范围和浓度分布, 在空气质量数值模拟和大气污染扩散问题的研究中倍受重视 (张强, 2001, 2003; 姜金华和彭新东, 2002; 何建军等, 2013; Reen et al, 2014)。因此, 大气边界层高度的研究对气象模式及空气污染数值模拟等有着重要意义, 是大气科学领域很多研究者关注的重点内容之一。

影响大气边界层高度的因素有很多, 边界层高度的发展不仅与大气层结稳定度等有关, 还受到下垫面的影响 (赵世强等, 2012; 涂静等, 2012)。目前大气边界高度的确定方法主要有三类 (赵鸣等, 1991; Liu and Ohtaki, 1997; 杜萍等, 2001; Hyun et al, 2005; 王耀庭等, 2012): (1) 基于廓线测量的方法, 主要有位温廓线法、风速廓线法、激光雷达资料反演法等; (2) 基于边界层基础理论的公式法, 主要有风速极值法、湍流能量法等; (3) 基于参数化和简单模型的数值模式法。其中, 风速极值法和湍流能量法在推导公式的过程中引入了很多的假设条件, 而这些假设在实际中很难满足, 导致其计算得到的边界层高度存在着较大的误差; 位温廓线法利用实际观测资料来确定大气边界层高度, 观测资料较难获得, 但其计算过程中没有引入假设, 计算结果比较可靠 (杜萍等, 2001; 张强等, 2004)。在风、温廓线数据充分的条件下, 常用位温廓线法来确定边界层高度 (张霭琛和田毅, 1990; Hägeli et al, 2000; 张鑫等, 2006); 参数化方法实用性较强, 但边界层高度变化的复杂性和模式各参数化方案的准确性还需进一步验证和改进 (王欣等, 2003)。

随着数值模式的广泛应用, 研究者试图用数值模拟的方法来表征大气边界层结构和湍流特征。WRF模式适用的尺度范围非常广泛, 是科学研究和业务预报中应用最为广泛的中尺度数值模式之一。近年来开展的很多科研工作进一步验证了WRF模式对边界层气象要素的模拟能力。安兴琴和吕世华 (2007)利用系留探空资料结合数值模式, 对比分析了兰州市冬季大气边界层温度场特征, 表明数值模拟结果能够反映温度场的基本特征; 缪国军等 (2007)利用WRF对兰州冬季风场和温度场进行了模拟, 结果表明WRF模式具有良好的稳定性; 王颖等 (2010)利用WRF模式不同边界层参数化方案模拟了兰州市冬季地面温度和风速的变化, 并将模拟结果与同期实测资料进行了对比分析, 验证了WRF模式模拟边界层气象场特征的能力; 徐慧燕等 (2013)利用WRF模式中不同边界层参数化方案对发生在长江中下游地区的3例暴雨进行了模拟实验, 并对降水总量和关键基本气象要素场的模拟结果与实测结果进行了统计检验; 张小培和银燕 (2013)对WRF模式四种参数化方案模拟复杂地形地区大气边界层中气象要素进行了评估; 王腾蛟等 (2013)利用WRF模式YSU和MYJ方案模拟了黄土高原丘陵地形条件下的夏季大气边界层结构, 结果表明YSU方案可更好地模拟出其对流边界层结构特征。利用WRF模式探讨大气边界层高度的研究也有很多, 陈炯和王建捷 (2006)利用WRF模式对北京地区夏季边界层日变化做了模拟研究, 认为高分辨率的模拟基本上可以反映出边界层的日变化特征; 张碧辉等 (2012)利用WRF两种边界层方案 (MYJ, YSU) 对沈阳冬季大气温度层结、低层风场、边界层高度等对污染物扩散有重要影响的气象要素进行了模拟; 何建军等 (2012)利用WRF模式Noah陆面过程参数化方案, 模拟了兰州地区冬季气象场, 并对模拟结果进行了检验; 马敏劲 (2010)采用2000年5月底至6月初在我国西北敦煌地区开展的“敦煌观测试验”所得到的探空资料和地面观测资料, 研究分析了试验期间观测到的超常规大气边界层高度, 并首次利用WRF模式对该极端大气边界层高度进行了模拟试验, 得到了与观测事实相符合的模拟结果, 验证了WRF模式对极端边界层高度的模拟能力; 赵世强等 (2012)利用激光雷达资料并结合WRF模式YSU和MYJ两种参数化方案估算兰州远郊榆中地区夏季边界层高度, 表明激光雷达反演边界层高度与WRF模拟边界层高度结果基本一致。国内西北地区大气边界层高度的研究主要集中在西北干旱区 (张强和王胜, 2008; 乔娟, 2009; 马敏劲, 2010; 韦志刚等, 2010; 何建军等, 2012; 赵建华等, 2013), 对于地处半干旱区的兰州地区, 大气边界层的研究主要集中在风温场 (安兴琴和吕世华, 2007; 缪国军等, 2007; 王颖等, 2010), 而针对冬季边界层高度的研究甚少。

本文利用WRF模式中三种不同参数化方案 (YSU、MYJ和ACM2) 对兰州市2005年1月2031日冬季大气边界层高度进行模拟, 分析典型背景条件下兰州市冬季边界层高度的日变化特征, 并将模拟结果与同期系留探空资料位温廓线法计算得出的边界层高度进行对比分析, 以探讨不同边界层参数化方案对边界层高度模拟结果的差异。

2 WRF模式设计方案 2.1 模拟期间天气形势

从2005年1月25-29日08:00(北京时, 下同)500 hPa环流形势 (图 1) 并结合实际观测资料可知, 20-28日期间兰州地区多处于平直的西风气流中, 受西北气流控制, 天气以晴为主。28日新疆上空形成了一个低温闭合中心, 该中心的温度槽与低压槽相配合, 形成了一次影响兰州的强降温降水过程, 30日以后气温开始缓慢回升 (图略)。

2.2 WRF模式模拟区域

利用WRFv3.4模式, 采用NCEP (1°×1°) FNL分析资料, 模拟区域为35.826 °N-36.3156 °N, 103.267 °E-104.367 °E, 区域中心点位于 (103.82 °E, 36.05 °N), 模式水平方向采用三重嵌套网格, 第一重网格区域为77×53, 格距25 km; 第二重网格区域131×91, 格距5 km; 第三重网格区域85×44, 格距1 km。垂直方向分为不等距的35层。图 2为第三重模拟区域地形高度, 其中西固二水厂位于36.12 °N, 103.61 °E。

2.3 模式参数设置

微物理过程参数化选取Lin et al (1983)的方案; 大气辐射参数化方案选取RRTM长波辐射方案 (1997)和Dudhia短波辐射方案 (1989); 积云对流参数化方案在第一、二重区域选取Kain-Fritsch (new Eta) 方案 (Kain, 2004), 第三重模式由于分辨率较高 (1 km) 而没有启动积云对流参数化方案; 陆面过程参数化采用热扩散方案 (1994)。模拟过程中微物理过程、大气辐射、积云对流及陆面过程等物理过程参数化方案保持一致, 但边界层参数化方案分别采用YSU、MYJ和ACM2方案 (Janjic et al, 1994; Hong et al, 2006; Pleim et al, 2007b), 对应近地面层方案分别选用Monin-Obukhov为基础的MM5相似方案、MYJ Monin-Obukhov方案和Pleim-Xiu surface layer方案 (Wang et al, 2010), 以探讨不同方案模拟边界层高度特征的差异。WRF模式参数配置如表 1所列。

2.4 边界层参数化方案介绍

WRFv3.4模式中提供了多种边界层参数化方案, 不同的边界层参数化方案的闭合方法和湍流交换机制存在很大的差异 (Mellor and Yamada, 1982; Hong et al, 2006; Pleim et al, 2007a; Pleim et al, 2007b; Wang et al, 2010; Hu et al, 2010; Xie et al, 2012)。本次模拟根据边界层不同的闭合方法选用YSU、MYJ和ACM2三种方案, 各个方案的闭合方法、湍流交换过程和边界层高度的确定 (王颖等, 2010) 详见表 2。

2.4.1 YSU方案

YSU方案 (Hong et al, 2006; Wang et al, 2010; 张碧辉等, 2012; Xie et al, 2012; 张小培和银燕, 2013) 是一阶非局地闭合的K理论 (梯度输送理论) 方案, 它是由MRF方案改进而来的, YSU方案中反梯度输送 (非局地通量) 项的量级比MRF方案中小, 使大气层结接近中性, 解决了MRF方案中反梯度输送项过大而导致的大气层结过于稳定的问题; YSU方案单独处理大气中的夹卷过程, 而不是把它作为边界层内部的混合过程, 从而增加了边界层中热力作用诱导的湍流混合, 降低机械 (动力) 强迫的湍流混合, 解决了在强风存在的条件下混合层中出现过度混合的问题; 同时, YSU方案还纠正了Hong et al (2006)的案例中PBL增长过快的问题。

YSU方案定义边界层高度h为逆温层中湍流通量最低值所在的高度, 即:

$ h=Ri{{b}_{\text{cr}}}\frac{{{\theta }_{va}}{{\left| U\left( h \right) \right|}^{2}}}{g[{{\theta }_{v}}\left( h \right)-{{\theta }_{s}}]}, $

(1)

式中: Rib_cr为临界里查逊数, U(h) 为h处的水平风速, θ_va为模式最底层的虚位温, θ_v(h) 为h处的虚位温, θ_s为近地面虚位温。

$ {{\theta }_{s}}={{\theta }_{va}}+{{\theta }_{T}}\left[ =b\frac{\overline{{{\left( w\prime \theta \prime v \right)}_{o}}}}{{{w}_{s}}} \right],\text{ } $

(2)

式中: θ_T为位温增量, w_s=u_*ϕ_m^-1为混合层速度 (u_*为表面摩擦速度), $\overline{{{\left( \text{w}'\theta {{'}_{v}} \right)}_{o}}}$为地表热通量, b=7.8。

首先, 不考虑θ_T(θ_s=θ_va), 根据公式 (1) 估算出h, 利用估算的h计算w_s和θ_T, 然后再根据公式 (1) 可求出边界层高度。

WRFv3.4模式YSU方案中, 不稳定层结时取临界里查逊数值为0;稳定层结时, 在陆地上, 将临界里查逊数值由0修正为0.25, 使得稳定层结下混合作用加强。

2.4.2 MYJ方案

MYJ方案 (Wang et al, 2010; Xie et al, 2012) 为Mellor-Yamada的2.5阶湍流闭合模型, 引入了湍流动能方程。MYJ方案边界层高度由湍流动能廓线决定, 定义边界层高度为湍流强度下降到临界值0.001 m²·s^-2时的高度; 该方案适用于所有稳定条件和弱不稳定条件的边界层, 但在对流边界层中误差较大; 由于MYJ方案考虑的物理过程比较复杂, 当不能满足其假设条件时, 模拟结果存在一定偏差。

2.4.3 ACM2方案

ACM2方案 (Xie et al, 2012; Pleim et al, 2007a, 2007b) 在ACM1方案的向上的湍流输送过程的基础上增加了局地输送部分, 它可以模拟由浮力作用引起的气块向上输送, 也可以模拟局地湍流交换。ACM2方案采用理查逊数方法计算边界层高度。

在稳定条件下, ACM2方案采用与YSU方案类似的方法诊断边界层高度, 即:

$ h=R{{i}_{\text{crit}}}\frac{\overline{{{\theta }_{v}}}U{{\left( h \right)}^{2}}}{g\left[ {{\theta }_{v}}\left( h \right)-{{\theta }_{v}}\left( {{z}_{1}} \right) \right]}, $

(3)

其中: θ_v为虚位温, z₁为模式最底层的高度, $\overline{{{\theta }_{v}}}$为在模式第一层和h之间的垂直平均虚位温, Ri_crit为临界体积里查逊数, Ri_crit=0.25时的高度定义为边界层高度。

不稳定条件下, 边界层高度的确定考虑了夹卷层的热力穿透和风切变 (不稳定情况下自由对流层顶部为浮力中性的夹卷层, 该层虽为稳定层结, 但由于热力穿透和风切变, 这里仍然有湍流)。Pleim (2007a, 2007b) 认为不稳定层结下PBL包括位于自由对流层以上的夹卷层, 边界层高度则为自由对流层和夹卷层之和。夹卷层的厚度定义为自由对流层顶与体积里查逊数达到一个临界值时所对应的高度之间的差值。在不稳定条件下, 先确定底部不稳定层结的高度 (Z_mix)。当θ_ν(z_mix)=θ_s时的高度即为对流不稳定层的高度,

$ {{\theta }_{s}}={{\theta }_{\nu }}\left( {{z}_{1}} \right)+b\frac{\overline{{{\left( w\prime \theta {{\prime }_{v}} \right)}_{o}}}}{{{w}_{m}}}, $

(4)

式中: b=8.5, $\overline{{{\left( \text{w}'\theta {{'}_{v}} \right)}_{o}}}$为地表运动热通量。然后根据理查逊数计算混合层高度, 即:

$ R{{i}_{b}}=\frac{g\left[ {{\theta }_{v}}\left( h \right)-{{\theta }_{s}} \right]\left( h-{{Z}_{\text{mix}}} \right)}{\overline{{{\theta }_{v}}}{{\left[ U\left( h \right)-U\left( {{Z}_{\text{mix}}} \right) \right]}^{2}}}, $

(5)

定义Ri_b=Ri_crit=0.25时所对应的高度为边界层高度。

3 系留探空资料选取 3.1 资料来源

兰州大学大气科学学院于2005年1月15日至2月2日期间在兰州西固区环境影响后评估研究中, 对西固区的污染气象和边界层特征进行了较全面的观测。大气边界层观测实验地点位于西固区二水厂, 本次研究采用实验期间系留探空资料 (AIR, TS-2A, 美国):观测要素为温度、气压、风向和风速, 观测时次为每天8次 (02:00, 05:00, 08:00, 11:00, 14:00, 17:00, 20:00, 23:00, 观测最大高度为400~800 m; 观测精度:温度0.2 ℃、气压1.0 hPa、风向5.0 °及风速0.25 m·s^-1。

3.2 实测资料计算边界层高度

利用试验期间系留探空资料, 采用位温廓线法, 根据边界层内位温廓线特征, 并结合前后时刻逆温发展的强度, 确定边界层高度。为了说明位温廓线法确定的边界层高度的合理性, 参考Dai et al (2011)在典型稳定和对流边界层下温度和位温梯度廓线特征 (由于系留探空资料位温观测是不连续的, 采用差分${{\left. \frac{\Delta \theta }{\Delta z} \right|}_{i}}=\frac{{{\theta }_{i}}-{{\theta }_{i-1}}}{{{z}_{i}}-{{z}_{i-1}}}$(i=2, 3, …, 21) 取代$\frac{\partial \theta }{\partial z}$作为位温梯度):稳定边界层以上$\frac{\partial \theta }{\partial z}$急剧减小; 在对流边界层以下温度近为常量, 在边界层顶温度随高度开始变化, 并且$\frac{\partial \theta }{\partial z}$明显增大, 图 3给出了不同层结条件下两种方法确定的边界层高度, 两者的结果基本一致。

由图 4可以看出, 边界层高度日变化的基本特征为白天太阳辐射增强, 地面吸收太阳辐射增温快, 夜间辐射冷却形成的逆温层逐渐受到破坏, 大气层结由夜间的稳定层结向不稳定发展, 湍流运动增强, 促进了大气边界层的发展, 边界层高度最高值出现在14:00左右, 可以达到600 m, 午后, 太阳辐射减弱, 大气边界层高度逐渐降低; 夜间, 地面辐射冷却降温, 形成逆温层, 大气层结稳定, 不利于湍流运动, 从而抑制了边界层向上发展, 逆温层高度在02:00左右最低 (约为80 m), 08:00以后, 由于太阳辐射增强, 地气之间温度梯度增大, 湍流增强, 边界层得到发展。本文利用观测资料位温廓线法结合位温梯度确定的边界层高度的日变化 (图 4b) 与杜萍等 (2001)计算的兰州大学和皋兰县冬季的边界层高度的日变化特征相似。

4 WRF模拟结果分析

模式初始条件和地形条件一致的情况下, 由于WRF模式不同边界层参数化方案采用的闭合方法、湍流交换机制及混合层高度的确定等存在很大的差异, 各参数化方案模拟的大气边界层高度存在很大差异。本文在对WRF模式不同参数化方案模拟的边界层高度的差异进行探讨的基础上, 通过统计分析给出模拟效果最优的边界层参数化方案。

4.1 兰州冬季PBLH的日变化特征

选取1月2327日期间WRF模式模拟的边界层高度进行分析, 在此期间没有大的天气过程, 可以代表兰州市冬季典型背景条件下边界层高度的变化特征。

从试验期间WRF模式三种边界层参数化方案模拟的边界层高度的日变化 (图 5) 可看出, WRF模式YSU、MYJ和ACM2三种边界层参数化方案模拟的PBLH的总体变化趋势相对来说是比较一致的, 均模拟出了兰州市冬季大气边界层高度日变化的基本特征。20:00至次日08:00, 大气层结处于稳定状态, 受兰州市山谷地形的影响, 逆温较强, 近地面湍流得不到充分发展, 边界层高度比较低, 夜间稳定边界层高度为150 m左右; 日出后, 太阳辐射增强, 湍流得到发展, 边界层高度逐渐抬升, 15:00左右达到最高, 持续一个多小时, 16:00以后湍流开始减弱, 边界层高度逐渐降低, 边界层趋于稳定。

由于不同边界层参数化方案所采用的闭合方法、湍流交换机制以及其确定PBLH的考虑不同, 各个参数化方案模拟所得的PBLH存在一定的差异: (1) 总的来说, 无论是夜间还是白天, MYJ方案模拟的边界层高度比YSU和ACM2方案均高, 这可能与MYJ方案对边界层高度的定义有关, MYJ方案定义湍流强度下降到临界值0.001 m²·s^-2时的高度为边界层高度, 而在边界层以上的夹卷层或残留层仍有零星湍流, 其强度可能超过临界值; (2) YSU和ACM2方案模拟的边界层高度的日变化特征比较相似, 不同的是: ACM2方案日间最大混合层高度高于YSU方案, 而夜间稳定边界层高度低于YSU方案。09:00-17:00, ACM2方案模拟的边界层高度高于YSU方案主要原因可能是ACM2方案在不稳定条件下边界层高度的确定考虑了夹卷层的热力穿透和风切变, 即在不稳定条件下, ACM2方案中边界层高度包含了自由对流层以上的夹卷层, 而YSU方案认为边界层不包括夹卷层; 17:00至次日09:00, 大气处于稳定层结, 由前述边界层参数化方案可知, 稳定层结时YSU和ACM2方案确定边界层高度的方法是类似的, 但在YSU方案的计算公式中考虑了θ_T, 增大了其计算的边界层高度 (Wang et al, 2010), 这是导致稳定层结时YSU方案模拟的边界层高度比ACM2的结果偏高的可能原因。

4.2 PBLH日变化特征对比分析

将WRF模式模拟的结果与同期系留探空资料位温廓线法确定的PBLH比较分析 (图 6), 以探讨WRF不同边界层参数化方案对兰州地区冬季大气边界层高度的模拟效果。

由图 6可以看出, 白天, 各参数化方案模拟结果整体较实测资料计算结果偏高, 尤其是在14:00左右, MYJ方案模拟的PBLH最高, ACM2方案次之, YSU方案最低; 模式模拟和实测资料计算的PBLH最高值均出现在14:00左右; 夜间 (23:00至次日08:00), MYJ方案和YSU方案模拟的结果均比实测资料计算值偏高, 而ACM2方案较低, 其中, MYJ方案模拟的PBLH的变化趋势与实测资料比较吻合, 但数值较实测值偏高; YSU方案的模拟结果比实测值略高, ACM2方案的结果较实测值略低。

为了定量分析三种边界层参数化方案的模拟效果, 下面将模式模拟的边界层高度与同期系留探空资料计算的边界层高度做统计检验 (选取2005年1月23-27日共5天40个样本)。选取平均绝对误差、平均相对误差、均方根误差和相关系数四个常用的统计量对模式模拟的结果进行检验 (三种边界层参数化方案均通过了α=0.01的显著性检验)。下面给出这4个统计量的计算方法:

$ MAE=\frac{1}{N}\sum\limits_{i=1}^{N}{\left| {{P}_{i}}-{{O}_{i}} \right|}, $

(6)

$ MAPE=\frac{1}{N}\sum\limits_{i=1}^{N}{\left( \left| {{P}_{i}}-{{O}_{i}} \right|/\left| {{O}_{i}} \right| \right)}, $

(7)

$ RMSE=\frac{1}{N}\sum\limits_{i=1}^{N}{{{\left( {{P}_{i}}-{{O}_{i}} \right)}^{2}}}, $

(8)

$ R=\frac{\sum\limits_{i=1}^{N}{\left( {{P}_{i}}-\bar{P} \right)-\left( {{O}_{i}}-\bar{O} \right)}}{\sqrt{\sum\limits_{i=1}^{N}{{{\left( {{P}_{i}}-\bar{P} \right)}^{2}}}}\times \sqrt{\sum\limits_{i=1}^{N}{{{\left( {{O}_{i}}-\bar{O} \right)}^{2}}}}}, $

(9)

式中: P_i为第i时刻的模拟值, O_i为第i时刻的观测值, P为模拟的平均值, O为观测的平均值, N为数据时间序列的长度。

表 3给出了不同参数化方案模拟的PBLH与同期西固系留探空资料用位温廓线法确定的PBLH的误差统计。

由表 3可知, ACM2方案模拟的兰州市冬季边界层高度与实测资料位温廓线法求得的PBLH之间的相关系数达到了0.91, 平均相对误差为49%, 优于YSU和MYJ方案。

综上所述, ACM2方案能较好地模拟兰州市冬季大气边界层高度。

5 结论

采用NCEP (1°×1°) FNL资料, 利用WRFv3.4模式模拟了2005年1月20-31日兰州市边界层高度。根据模拟结果, 比较分析了YSU、MYJ和ACM2三种边界层参数化方案模拟的边界层高度, 并将WRF模式模拟结果与同期系留探空资料位温廓线法计算的边界层高度进行对比分析。得到如下主要结论:

(1) WRF模式三种不同的参数化方案 (YSU、MYJ和ACM2) 模拟的边界层高度总体变化趋势一致, 均模拟出了兰州市冬季PBLH日变化的基本特征。

(2) 由于不同边界层参数化方案所采用的闭合方法和混合层高度的确定方法等不同, 各方案模拟的边界层高度存在一定的差异:总体来说, 无论是夜间还是白天, MYJ方案模拟的边界层高度明显偏高; YSU和ACM2方案模拟的边界层高度日变化特征比较相似, 两者主要区别在于: ACM2方案日间最大混合层高度高于YSU方案, 而夜间稳定边界层高度低于YSU方案。

(3) 统计分析结果显示, ACM2方案能较好地模拟兰州市冬季大气边界层高度。