1 引言

青藏高原平均海拔4000 m以上(张镱锂等, 2002)，是中国最大、世界上海拔最高的高原。其特殊的动力和热力作用不仅会影响局地地-气相互作用，还对大气环流、季风和气候变化等产生重要影响(吴国雄等，2004；徐祥德等, 2006；朱玉祥等，2007；马耀明等, 2014)。青藏高原的动力和热力作用通过大气边界层传递到自由大气。因此，对青藏高原大气边界层特征的分析显得尤为重要。

青藏高原具有独特的大气边界层特征(徐祥德等, 2001; 马耀明等, 2002)。青藏高原大气边界层内风速随高度升高呈现多极值分布(李家伦等, 2000)。干季，青藏高原大气边界层往往出现逆湿现象(Zuo, et al, 2005)，湿季，则逆湿现象消失。由于海拔高、日照强，在青藏高原往往可以观测到极大的太阳总辐射和净辐射(徐祥德等, 2006)。这也导致青藏高原上湍流交换强，大气边界层发展旺盛(周文等, 2018)。李茂善等(2011)研究发现，青藏高原那曲地区干季边界层高度为2211—4430 m, 雨季边界层高度为1006—2212 m，明显高于平原地区的1000—1500 m (赵鸣等, 1992)。

目前，针对青藏高原大气边界层研究已经做了许多的工作(王介民, 1999; Wang, et al, 2016)。如第一次青藏高原大气科学考察试验(QXPMEX)(翁笃鸣等, 1984)、第二次青藏高原大气科学考察试验(TIPEX)(徐祥德等, 2006)、全球协调加强观测计划亚澳季风之青藏高原试验研究(马耀明等, 2006)，以及正在进行的第三次青藏高原大气科学试验等。但由于恶劣的观测环境和高昂的观测费用等原因，青藏高原大气边界层研究多局限于近地层气象要素和湍流通量的分析，对气象要素的垂直变化特征研究较少。

随着数值模式的发展和计算机计算能力的提高，数值模式成为大气科学研究的重要工具。数值模式中采用边界层参数化方案来表征大气边界层的变化特征。Krogsæter等(2015)采用WRF模式中5种边界层参数化方案对北海离岸风进行模拟, 发现MYJ方案对海洋边界层的模拟与观测更为一致。Gómez-Navarro等(2015)对不同边界层参数化方案的风场模拟效果进行了对比分析，发现非局地边界层方案能够更好地解析复杂地形下边界层风的特征。Banks等(2016)采用WRF模式中8种边界层参数化方案对不同天气形势下雅典城市下垫面近地层气象要素和边界层高度进行研究，发现非局地方案能够更准确地模拟边界层高度。中国也开展了相关研究。张碧辉等(2012)研究发现，YSU方案边界层混合作用较强，模拟边界层高度较高，有利于污染物的垂直扩散。张小培等(2013)对复杂地形下大气边界层特征进行模拟分析，发现ACM2方案能更好地模拟边界层气温、风场变化特征。王成刚等(2017)也发现ACM2方案对晴天和阴天条件下温度和比湿的模拟效果更好。但是，由于观测资料的缺乏，不同边界层参数化方案在青藏高原的适用性研究较少(李斐等, 2017)。本研究基于第三次青藏高原大气科学试验边界层观测结果，对青藏高原那曲地区边界层参数化方案的适用性进行研究。

那曲地区位于青藏高原中部，平均海拔在4500 m以上，主要植被类型为高寒草甸、高寒草地和灌丛(杨秀海等, 2008)，是青藏高原的典型下垫面。文中利用WRF模式4种不同参数化方案(YSU、MYJ、ACM2和BouLac)对青藏高原那曲地区2015年5月30日—6月30日大气边界层特征进行模拟，分析那曲地区边界层气象要素的日变化特征，并将模拟边界层风、位温、相对湿度垂直廓线与探空资料观测结果进行对比，分析WRF模式不同边界层参数化方案对青藏高原边界层特征模拟效果的影响。边界层参数化方案在不同下垫面、地形和天气条件下的模拟效果存在一定的不确定性，通过更具有代表性的个例来评估不同参数化方案的模拟效果，对于边界层参数化方案的改进和发展具有一定的现实意义。

2 模式和数据 2.1 模式设置

采用中尺度天气预报模式WRF对青藏高原那曲地区边界层特征进行数值模拟。模拟采用4层双向嵌套，从外层向内层的格点数分别为145×160、289×169、229×142和550×523，分辨率依次为27、9、3和1 km(图 1)。模式输入资料为ERA-Interim 0.25°×0.25°再分析资料，逐6 h输入一次，地形资料为MODIS(30″)格点资料。模式参数方案设置如表 1所示。

WRFv3.7.1模式提供多种边界层参数化方案。不同方案模拟的湍流混合强度不同，导致边界层内温度、风向风速模拟结果存在明显差异(Xie, et al, 2012)。依据不同闭合方法和湍流混合机制，文中采用两种局地方案(MYJ、BouLac)和两种非局地方案(YSU、ACM2)分别进行模拟，不同边界层参数化方案闭合方法、湍流混合及边界层高度定义如下：

YSU方案(Hong, et al, 2006)是由MRF方案改进的一阶非局地闭合的K理论方案。其湍流扩散方程为

(1)

式中，C为诊断变量，t为时间，z为高度，K_c为湍流扩散系数，γ_c为局地梯度修正项(反梯度项)，为逆温层湍流通量，h为边界层高度。YSU方案中用反梯度输送(非局地通量)项表示由非局地梯度引起的湍流交换，并减小反梯度项量级，以此解决MRF方案中大气层结过于稳定的问题。YSU方案定义边界层高度为逆温层中湍流通量极小值所在层的高度，稳定层结时取临界理查森数值为0.25的高度，不稳定层结时取临界理查森数值为0的高度。

MYJ方案(Janjić, 2002)为Mellor-Yamada 2.5阶局地湍流动能方案。其湍流动能方程为

(2)

式中，e为湍流动能，ρ为空气密度，u为纬向风分量，v为经向风分量，w为垂直速度，z为高度，β为浮力系数，θ为位温，ε为耗散项，带“′”量为相应物理量的脉动量。MYJ方案通过计算湍流动能来确定湍流扩散系数，适用于稳定和弱不稳定大气边界层(张碧辉等, 2012)。由于MYJ方案考虑的物理过程较为复杂，当不满足其假设条件时，模拟结果存在一定偏差。MYJ方案定义湍流强度下降到临界值0.2 m²/s²时的高度为边界层高度。

ACM2方案(Pleim, 2007)是一阶闭合方案，其是在ACM1方案的基础上增加了局地湍流输送，可以模拟由浮力作用引起的非局地向上输送和局地的湍流交换。其湍流扩散方程为

(3)

式中，C_i为标量C在模式第i层的混合率，f_conv为调控局地混合和非局地混合的比例系数，M_u为非局地向上对流混合率，M_di为局地向下混合率，Δz_i为模式第i层厚度。稳定或中性层结时，ACM2关闭非局地而采用局地湍流输送。边界层内湍流交换系数由局地理查森数湍流交换公式和相似理论计算得到，不稳定层结中，湍流交换系数取两者间较大值。ACM2方案采用理查森数方法计算边界层高度。稳定层结时取临界理查森数值为0.25的高度；不稳定层结时边界层高度则为自由对流层和夹卷层高度之和。

BouLac方案(Bougeault, et al, 1989)为1.5阶局地湍流动能方案，其湍流动能方程同式(2)。与MYJ方案不同，BouLac方案在湍流交换系数和混合长的计算中考虑了地形引起的湍流混合。对陡峭地形下，该方案通过湍流动能廓线来确定边界层高度，边界层高度定义为湍流动能减至0.005 m²/s²时的高度。

2.2 数据来源

验证WRF模拟结果所用观测资料来源于“第三次青藏高原大气科学试验”在青藏高原那曲地区的5个观测点，各观测点地理位置及海拔高度如表 2所示。观测点平均海拔高度4500 m以上，地表植被类型为高寒草甸。观测区域地形平坦，所测资料基本上可以代表青藏高原那曲地区边界层特征。各观测点主要观测基本气象要素、辐射以及湍流通量等。观测场内建有2—3 m高的边界层观测塔，上面搭建了温、湿度传感器(VAISALA HMP45C-GM)和风向、风速传感器(METONE 034A/B)、三维超声风速温度仪(CSAT3)和CO₂/H₂O红外气体分析仪(Li-Cor, LI-7500)。三维超声风速温度仪用于测量u 、v、w3个方向风速脉动和声学虚温，采样频率10 Hz。CO₂/H₂O红外气体分析仪用于测量大气中CO₂和水汽的密度，采样频率10 Hz。并采用热流量板测量土壤热通量(HFP01)。涡动观测资料采用EddyPro 4.0软件进行处理并进行过质量控制。此外，采用那曲观测点探空资料对垂直方向风向风速、位温廓线的模拟进行验证。模拟结果的定量评估采用平均偏差(MB)、相关系数(R)和均方根误差(RMSE)。

3 结果分析 3.1 模式结果验证

利用2015年6月青藏高原那曲地区5个观测点(安多、班戈、比如、嘉黎和那曲)的常规气象观测资料，与4种边界层参数化方案模拟的近地面气温、地表温度和风速进行对比分析，以验证模式结果，并分析青藏高原那曲地区的气象要素特征。YSU、MYJ、ACM2和BouLac边界层参数化方案模拟的2 m气温平均日变化与观测值的对比如图 2所示。各方案模式性能统计分析见表 3。安多、班戈、比如、嘉黎和那曲站2015年6月平均气温分别为8.05、8.65、9.78、7.76和10.04℃。观测气温07时(北京时，下同)达到最低值，17时达到峰值。各方案模拟的气温与观测变化趋势一致，谷值和峰值出现时间相同。但在安多、班戈、嘉黎和那曲站，各方案模拟的气温偏低。这可能是因为模式地形高于实际地形。YSU方案增加边界层中热力强迫引导的湍流混合，近地层气温最低。MYJ方案模拟的气温更接近观测值，与观测值的相关系数为0.96。与Xie等(2012)研究发现的局地方案模拟气温负偏差较大不同，相比于非局地方案(YSU和ACM2)，局地方案(MYJ和BouLac)在那曲地区模拟的气温负偏差较小。这可能是因为相较于城市下垫面，那曲地区下垫面类型较为均一，局地湍流动能方案反而更适用。

各方案模拟的地表温度明显低于观测。尤其是安多站，正午地表温度负偏差达10℃左右。Jiang等(2005)的研究同样发现青藏高原地表温度模拟负偏差的问题。Zhuo等(2016)的研究指出，地表与大气的感热输送偏差可能是导致地表温度负偏差的原因。相对于其他方案，BouLac方案在湍流扩散系数的计算中考虑了地形影响产生的湍流交换，模拟的地表温度负偏差较小。比如站，BouLac方案模拟的地表温度平均偏低0.41℃，与观测值的相关系数为0.96。

3.2 能量交换特征

陆-气能量交换对边界层气象要素的变化有很大影响。能量平衡各分量月均日变化特征如图 4所示。各方案模拟值的平均偏差、相关系数和均方根误差统计见表 5。安多、班戈、比如、嘉黎和那曲站能量平衡闭合率分别为0.92、0.87、0.82、0.80和0.90。Wilson等(2000)的研究发现，涡动相关系统观测的能量平衡闭合率为0.53—0.99。导致能量不闭合的原因主要有：没有考虑复杂地形下平流的作用、涡动观测对感热和潜热通量的低估、通过热流板观测热通量间接计算土壤热通量的估算误差等(Feng, et al, 2016)。

能量平衡各分量在日出后增大，午后达到最大值。净辐射14时达到峰值650 W/m²左右。安多、班戈和比如站，净辐射模拟值与观测一致。那曲站，净辐射模拟值低于观测。这是因为那曲站向下长波辐射模拟偏低。各方案模拟的净辐射差异不大。感热通量日变化范围为-30—360 W/m²。日出后陆地升温，陆-气温差逐渐增大，感热通量升高，08时许由负值转为正值。日落后地面辐射冷却降温，陆-气温差减小，感热通量降低。安多和班戈站感热通量观测值15时分别达到峰值353和251 W/m²。比如、嘉黎和那曲站感热通量观测值14时分别达到峰值174、142和274 W/m²。感热通量与陆-气温差存在较显著的正相关(Liu, et al, 2015)。安多站地表温度高于其他站，陆-气温差大，地面向大气输送的感热通量多。模式模拟感热通量偏大。MYJ方案在班戈站感热通量模拟正偏差达50 W/m²。这是因为复杂地形下模式地形过于平坦，次网格地形对环流的影响减小，风速偏大，陆气间热量交换偏大(Jiménez, et al, 2013)。潜热通量的日变化范围为0—280 W/m²。安多、班戈、比如、嘉黎和那曲站15时分别达到峰值118、61、276、200和188 W/m²。班戈站土壤含水量低，陆-气水汽压差小，潜热输送少。青藏高原绒布河谷春末夏初地表干燥，空气含水量低，仲雷等(2007)的研究中潜热通量峰值仅有54.47 W/m²。MYJ方案近地层湍流交换强，模拟潜热通量正偏差大。土壤热通量白天为正值，地面吸收热量升温，夜间为负值，地面释放热量加热大气。土壤热通量12时前后达到峰值，比净辐射、感热通量和潜热通量提前约2 h。白天，模式模拟土壤含水量低于观测，模拟土壤热通量偏低。相比其他方案，ACM2方案模拟的土壤热通量平均偏差较小，与观测值的相关系数最高。ACM2方案考虑了浮力作用引起的非局地向上输送和局地的湍流交换，能够更好地模拟云层状况(王成刚等，2017)，向下短波辐射和长波辐射的模拟偏差较小，地面接收热量更接近观测。

由于青藏高原海拔较高，吸收太阳辐射多，各观测点净辐射明显高于低海拔地区的550 W/m²左右(Xu, et al, 2016)。青藏高原下垫面的感热通量比沙漠(李建刚等, 2012)和湖泊(许鲁君等, 2014)大，潜热通量比两者小。有效能量(R_n-G)主要分配为感热通量。MYJ方案对感热通量和潜热通量模拟的正偏差较大，BouLac方案对感热通量模拟的正偏差较小。

3.3 边界层结构特征

利用那曲地区2015年6月1—10日的探空资料，验证4种边界层参数化方案模拟的边界层结构特征和气象要素垂直分布特征。探空资料的观测频率为每天2次，分别在08和20时进行观测。图 5为观测期间垂直风廓线与YSU、MYJ、ACM2和BouLac方案模拟风廓线的对比。观测的近地层主导风向为西风，距地面300 m左右高度后偏西北风，800 m左右高度风向转变为西南风。各方案模拟的近地层主导风向为西南风，ACM2方案模拟的近地层风向偏差约35°。Boulac方案模拟的近地层风向偏差较小，约为20°。随着高度的升高，模拟风向向西北风转变，地面高度500 m左右为西北风，900 m左右风向又转变为西南风。1500—3000 m模拟风向偏南，3000 m以上模拟风向与观测一致。模式模拟风速比观测大。BouLac方案模拟的近地层风速正偏差最小，其次为MYJ方案，YSU方案模拟的近地层风速正偏差最大。风速模拟的正偏差主要是经向风模拟偏大导致的。

图 6对比分析了08和20时位温廓线观测值和YSU、MYJ、ACM2和BouLac方案模拟值。观测期间，青藏高原那曲地区平均日出时间为07时前后，日落时间21时前后。观测的08时平均位温廓线仍表现为典型的稳定边界层层结。夜间，地面辐射冷却降温，近地面存在逆温层。各方案模拟出夜间的稳定边界层，但模拟的2000 m以下边界层偏冷。BouLac方案模拟的冷偏差较小，约为3.2 K，ACM2方案模拟的冷偏差较大，约为4.3 K。地面高度2700—4300 m，模式模拟位温高于观测值。夜间稳定边界层高度约为200 m，MYJ方案模拟的边界层高度与观测一致，其他方案模拟的边界层高度比观测高约100 m。

图 6b表现为对流边界层，近地层为超绝热逆温层。地面高度3500 m以下，模拟位温低于观测，最大负偏差为2500 m左右的5.2 K。地面以上3500 m左右为混合层顶部的卷夹层，各方案对卷夹层的模拟均存在偏差。李斐等(2017)的研究指出，当边界层顶部逆温层弱时，模式对边界层位温分布模拟效果较好。距地面高度2100 m以下，MYJ方案模拟位温冷偏差最小，其次为BouLac方案，ACM2方案模拟的冷偏差最大。距地面高度4200—5300 m，模拟位温高于观测值。随着高度升高，模拟偏差减小。白天对流边界层高度约为3400 m，各方案对边界层高度的模拟偏低。

模式对边界层的模拟偏冷，对4000 m左右大气的模拟偏高。夜间模拟的冷偏差小于白天，表明湍流交换较弱时，模式模拟效果较好。无论是白天还是夜间，局地方案的模拟效果均优于非局地方案。在下垫面类型相对均一的青藏高原那曲地区，基于湍流动能理论的边界层方案模拟效果比基于K理论的边界层方案模拟效果好。ACM2和YSU方案较强的非局地作用使得边界层气温模拟冷偏差较大。

边界层内的水汽分布会直接影响云的形成，并对地表能量吸收产生影响。08时，近地面相对湿度为27%，随着高度的升高，相对湿度减小(图 7)。距地面高度500—2300 m，相对湿度随高度的升高升高。距地面高度2300—3200 m，相对湿度随高度的升高基本不变，之后随着高度的升高降低。各方案对1500 m以下相对湿度的模拟偏高，BouLac方案对相对湿度模拟的正偏差较小，ACM2方案模拟的正偏差最大。WRF模式对近地层湿度模拟偏高，这可能是模式对青藏高原下垫面(如土壤含水量、植被覆盖率)描述失真导致的(Yang, et al, 2018)，陆-气较强的水汽输送也会导致边界层湿度模拟偏高。模式对1500 m以上相对湿度的模拟偏低。距地面1800—3500 m高度，ACM2方案模拟的相对湿度负偏差最小。20时，距地面3500 m高度以下，相对湿度随高度的升高降低。距地面3500—4500 m高度，相对湿度随高度的升高变化不大，之后随着高度的升高降低。距地面3000 m高度以下，模式对相对湿度的模拟偏高。MYJ和BouLac方案模拟的相对湿度正偏差较小，ACM2方案模拟的正偏差较大。YSU、MYJ、ACM2和BouLac方案模拟的相对湿度正偏差在2500 m左右高度达到最大值，分别为45 %、42 %、40 %和23 %。模式对陆-气的水汽交换模拟偏强。相比于其他方案，BouLac方案模拟的相对湿度与观测更为一致。ACM2方案对边界层内水汽交换的模拟偏差较大，这与李斐等(2017)的研究一致，与王成刚等(2017)的研究结果相反。

4 结论

利用“第三次青藏高原大气科学试验”在安多、班戈、比如、嘉黎和那曲站的观测资料，采用中尺度天气预报模式WRF中4种边界层参数化方案对青藏高原那曲地区边界层风、温、湿和垂直结构特征进行数值模拟，分析不同边界层参数化方案在该地区的适用性。得到的主要结论如下：

(1) WRF模式对青藏高原那曲地区近地面气温和地表温度的模拟存在冷偏差。相对于非局地方案(YSU和ACM2)过强的混合作用，局地方案(MYJ和BouLac)在那曲地区模拟的气温冷偏差较小。MYJ方案对近地层气温的模拟与观测最为一致，与观测值的相关系数为0.96。BouLac方案对地表温度模拟效果最好。

(2) 在青藏高原那曲地区，地表吸收的能量主要以感热的形式释放。模式模拟的感热通量和潜热通量偏大。复杂地形下，模式忽略了次网格地形对环流的影响，风速模拟偏大，陆-气热量交换偏大。BouLac方案考虑了地形引起的湍流交换，对感热通量模拟的正偏差较小。从能量平衡的角度看，温度模拟的负偏差主要来源于向下长波辐射模拟偏低以及感热通量和潜热通量交换过强，这与Chen等(2017)的研究结果一致。

(3) WRF模式模拟的边界层风垂直分布特征与观测一致，但模拟风速比观测大，主要是经向风模拟偏大导致的。BouLac方案模拟的近地层风速正偏差最小。在下垫面类型相对均一的那曲地区，基于湍流动能理论的边界层方案比基于K理论的边界层方案模拟的位温和相对湿度与观测更为一致。

对于不同气象要素，不同边界层参数化方案的模拟效果不同。总体来看，考虑了地形因素的BouLac方案对那曲地区温度、位温和水汽垂直分布的模拟效果较好。但本研究仅对2015年6月无降水天气的边界层特征进行了对比分析，对青藏高原边界层参数化方案的选取具有一定的参考意义，但对季风期降水时的边界层特征和地形复杂的其他地区尚有待进一步分析。

致谢: 感谢“第三次青藏高原大气科学试验”项目组提供的相关数据及某些数据的解释说明。尤其感谢安多、班戈、比如、嘉黎和那曲站相关科研人员提供的观测数据。