青藏高原夏季边界层再分析资料的偏差分析及订正

马敏劲，陈 玥，康国强，赵侦竹，黄万龙，谈昌蓉,丁 凡

(1.甘肃省干旱气候变化与减灾重点实验室，甘肃 兰州 730020；2.兰州大学大气科学学院，甘肃 兰州 730000；3.兰州理工大学计算机与通信学院，甘肃 兰州 730050)

引 言

大气边界层作为自由大气与下垫面之间互相作用的媒介，其内部的动力、热力过程对地气系统能量平衡具有重要意义[1]。一些重要的天气现象如降水、雾、霜等，特别是沙尘暴、暴雨等灾害性天气现象的形成与边界层过程密切相关。青藏高原平均海拔高度4 km以上，大地形抬升增强了地面对太阳辐射的吸收，高原热力作用显著，边界层内对流系统发展旺盛[2-3]。另一方面，在地表感热不是很强的冬季，高原上空观测到厚度为5 km左右的超厚对流边界层[4]，边界层顶部接近9 km，基本达到对流层顶的位置。高原对流发展旺盛直接影响天气系统，在西风带作用下传输影响至下游地区，甚至整个东亚[5-10]。

大气边界层高度是衡量对流发展的重要参量，其精准确定需要垂直加密观测，迄今为止在高原地区已经开展了数次大型观测试验，研究表明青藏高原边界层显著高于平原地区[11]，且季节变化明显，旱季(11月至次年3月)边界层在地面强感热通量作用下发展深厚，最高可达4500 m左右，雨季(5—9月)边界层发展由增强的潜热通量主导，最高在2500 m左右[12-13]。但由于观测环境恶劣，青藏高原常规气象观测站点稀疏，且站点无线电探空大部分限定在特定气压层内，垂直分辨率不足，导致高原大气边界层研究大多局限于近地面层通量和湍流特征，物理量垂直分布研究明显不足[14]，限制了高原超高的大气边界层发展机制及其对天气系统影响的研究[15]。

再分析资料具有覆盖性好、垂直分层高等特点，是研究青藏高原气象要素特征的重要资料，如MERRA是青藏高原及其周边地区较为理想的高分辨率再分析资料[16-17]，ERA-Interim数据可作为高原多年冻土区地表温度的有效代用资料[18]。也有研究指出在高原地区，再分析资料ERA-40和JRA-25的适用性要好于NCEP[19-21]，风场上ERA-Interim和CFSR资料总体好于JRA-55资料[22]，而ERA-Interim和MERRA分别能较好地刻画夏季水汽输送情况[23]。

高原大地形引起的山地复杂下垫面，通过陆-气相互作用形成高原上空复杂的大气涡动特征和热力结构变化[24]，虽然再分析资料一定程度上可以弥补高原地区观测资料缺乏的缺陷，但某些物理量场存在较大误差，有必要对再分析资料进行订正。如NCEP再分析资料在强对流环境分析中水汽参数与探空观测资料差异较大，大气边界层内水汽量需要观测资料进行误差订正[25]，山地风场需要利用比值法进行误差订正[26]。同样，再分析资料对物理过程刻画不足，需引入物理框架进行约束。在一定模式参数化方案的本地化条件下[27-28]，数值模式对于边界层结构的模拟能力，尤其极端大气边界层高度的数值模拟有较大改进[29-32]。

本文首先探析不同再分析资料在夏季青藏高原边界层的适用性，然后选择适用性好的再分析资料作为初始场开展数值模拟，对比分析不同参数化组合方案对高原大气边界层高度的模拟效果，利用模式中物理框架约束作用订正再分析资料在青藏高原边界层中的分析误差，以期深化再分析资料在青藏高原应用的认识。

1 数据与方法

1.1 数据来源

数据时段为2017年夏季(6—8月)，研究区域为青藏高原及其周边地区(65°E—110°E，25°N—43°N)。ERA-Interim由欧洲中期天气预报中心ECMWF (European Centre for Medium-Range Weather Forecasts)提供的新全球大气再分析资料[33]。ERA-Interim采用四维变分同化，并结合改进的湿度分析、卫星数据误差校正等技术，实现了再分析资料质量的提升和满足多种网格精度的需求。JRA-55是日本气象局(JMA)提供的1958年以后55 a的高质量同质性气候数据集，其时间跨度在第三代再分析资料中最长[34]。ERA-Interim和JRA-55的空间分辨率均为0.25°×0.25°，时间分辨率为6 h。卫星数据是NASA发布的高分辨率再分析资料MERRA-2逐小时单层同化诊断数据集(M2T1NXSLV)，空间分辨率为0.5°×0.625°，时间分辨率为1 h[35]。MERRA-2(the modern-era retrospective analysis for research and applications, version 2)是NASA基于第5版戈达德地球观测系统模型(GEOS-5, the goddard earth observing system model, version 5)及其大气同化系统(ADAS, atmospheric data assimilation system)进行的卫星大气再分析数据。由于长时间序列的探空观测资料有限，且垂直分辨率较低，为了减小观测数据误差带来的不确定性，选择NOAA和Wyoming大学提供的两套探空资料同时进行检验，并取二者的均值来检验结果。其中，来自NOAA的数据提供了高原上空6个测站资料，分别为西藏拉萨、那曲、昌都站，青海玉树站，四川甘孜站和甘肃合作站，Wyoming大学提供了除合作外的以上5个测站。图1为研究区域及6个观测站点示意图。

图1 青藏高原地形及Wyoming和NOAA资料的探空站点分布

文中附图涉及地图均基于国家测绘地理信息局标准地图服务网站下载的审图号为GS(2016)2885号的标准地图制作，底图无修改。

1.2 数值模拟参数化方案

对比的边界层参数化方案包括YSU、MYJ和ACM2方案，微物理方案包括Purdue Lin和WSM6方案，积云对流方案包括Kain-Fritsch(new Eta)和Betts-Miller-Janjic(简称“BMJ”)方案(表1)。

表1 试验中对比的模式参数化方案简介

1.3 模拟结果检验方法

模拟的气象要素结果主要利用平均偏差(Bias)、均方根误差(RMSE)和误差标准偏差(STDE)进行评估[49-50]，计算公式如下：

(1)

(2)

(3)

其中：Si表示模拟值；Oi代表观测值；n是样本量。在评估模式模拟效果时，通常会优先考虑STDE值，这种假设来源于一些事实，因为它代表了模式误差的离散程度。本质上，STDE值偏小代表误差接近于一个常数，即使此时Bias和RMSE值都比较大，但依然能说明模式的物理过程正确。相反，即使Bias和RMSE值都很小，但STDE是一个大值，则表示误差是随机分布的，说明模式缺少物理基础。

再分析资料通过加入模式物理框架约束后进行订正，因没有足够时长和较高垂直分辨率的探空观测资料来计算边界层高度，再分析资料的订正效果通过同高度的气象要素和边界层高度的季节分布来分析。

2 再分析资料在青藏高原边界层中的适用性

GUO Jianping等[51]研究中国北方及青藏高原的对流天气过程及其对下游天气的影响时，认为ERA-Interim再分析资料计算的大气边界层高度在中国北方最接近探空资料的计算结果，但在高原地区的适用性没有进一步探讨。以下分析ERA-Interim、JRA-55和NASA MERRA-2再分析资料确定大气边界层高度的关键气象要素，利用青藏高原的探空观测资料进行对比检验，探讨再分析资料在青藏高原地区边界层研究中的适用性。

2.1 探空检验

利用探空观测资料检验不同再分析资料在青藏高原不同高度上的气象要素，并用这些气象要素来确定大气边界层发展高度。虽然理查逊数法常用于模式参数化方案中，但高原上缺乏与温湿要素同一层次的风场资料及局地地表粗糙度参数化信息，使得边界层高度计算结果稳定性及可靠性较差。此外，探空站点集中于高原东南部，空间上不足以代表整个高原且垂直分辨率(100 hPa)低，无法捕捉到日间边界层顶的变化特征，夜间逆温位置计算不准确使稳定边界层高度存在严重高估。因此，首先选择使用边界层内特定等压面上的气象要素来讨论再分析资料在青藏高原地区的适用性。

由表2可知，除MERRA-2在高层缺失数据外，JRA-55和ERA-Interim都有完整数据，3套再分析资料都可以较好地反映青藏高原边界层内的温度(T)特征，而水平风场尤其是风向的描述效果较差。从偏差来看，3套再分析资料均低估青藏高原东南部边界层内露点温度(Td)，且除MERRA-2的500 hPa数据外，温度也存在普遍低估。需要注意的是，再分析资料的地面风向(DIR)、风速(SPD)偏差可能包含了探空资料误差，因为2组探空数据在近地面包含大量非缺测的0值。但是除近地面外，再分析资料的高空风速误差也远远大于其他变量，而探空资料在高层具有良好的连续性，再分析资料对于青藏高原边界层的风速描述不够准确。MERRA-2数据仅包含特定层次上的要素信息，所以可以认为该资料对高原边界层内部特征的描述能力有限。虽然总体来看JRA-55和ERA-Interim再分析资料在青藏高原东南部的适用性相当，但是JRA-55的露点温度偏差很大，且STDE值远大于另两套数据，说明其误差分布随机，即该数据不能很好地反映青藏高原边界层内的湿度场特征。因此夏季青藏高原东南部边界层适用性最高的再分析资料可以认为是ERA-Interim。

表2 再分析资料在青藏高原大气边界层气象要素的探空检验

再分析资料提供一日4个时次(00:00、06:00、12:00及18:00，世界时，下同)资料，探空观测提供一日2次(00:00和12:00)资料，数据不能完全支持分析大气边界层发展变化，缺失较多细节，中尺度数值模拟能为高原大气边界层发展的研究提供更多详细信息，而且模式参数化方案为再分析资料描述青藏高原大气边界层提供物理过程的约束，因此，以下借助模式模拟结果进一步分析青藏高原大气边界层特征。

2.2 ERA-Interim的边界层特征分析

分析较长时间边界层高度特征的最佳方法是由VOGELEZANG等[52]最初提出的临界理查逊数法，该方法对稳定边界层和对流边界层皆适用，且不强烈依赖于观测的垂直分辨率。理查逊数Ri是与浮力相关的湍流切变与机械切变的比值，其计算公式如下：

(4)

其中：z(m)为高度，z0(m)为地面海拔高度；θ(K)为位温，为了能同时描述日间对流边界层和夜间稳定边界层特征，用虚位温θv(K)代替θ；u，v分别为水平经向和纬向风速(m·s-1)；b是常数，u*是表面摩擦速度。

利用边界层气象要素检验整体结果最好的ERA-Interim再分析资料，得到其季节平均空间分布，如图2所示，高低值区分界线与青藏高原地形边界契合，在分界处高度梯度明显。夏季与冬季青藏高原及其北部的边界层高度分布呈现高低相反特征且高度梯度显著，夏季高原上空边界层高度在800～1000 m之间，高原北部边界层高度在1400 m以上；冬季高原上空边界层高度在1200～1600 m之间，高原北部边界层高度在400 m以下。高原上空边界层高度春季最高，高值中心出现在高原西侧，由西南向东北逐渐降低，与高原北部高度梯度较小。秋季高原上空与北部地区的边界层高度梯度最小，高原边界层高度略高于高原北部。从季节变化来看，高原上空边界层高度夏低冬高，同时春、夏季边界层高度日变化清晨高、傍晚低(图3)，与以往高原站点观测规律并不一致[11,13]。因此，具有较高空间分辨率的ERA-Interim再分析资料在青藏高原对边界层高度的刻画存在较大问题，以下将利用模式物理框架进行约束，观察边界层高度计算结果的改进。

图2 ERA-Interim再分析资料得到的边界层高度季节平均空间分布(单位：m)

图3 ERA-Interim再分析资料得到的季节平均边界层高度的日变化时间序列

3 模式物理框架的约束

3.1 试验设置

数值模拟选取探空资料站点集中的高原东南部地区开展试验，先探讨数值模式的物理过程参数化方案的适用性。WRF version 3.8中边界层参数化方案提供了13种方案，针对青藏高原的参数化方案选择成为数值模拟的一个重要问题，通过敏感性试验得到较好适用于本地的方案组合。

本次模拟水平方向采用三重网格嵌套(图4、表3)，第一重网格距为27 km，经向格点数163，纬向格点数88；第二重网格距为9 km，经向、纬向格点数分别为367、160；第三重网格距为3 km，经向、纬向格点数分别为367、226，第三重网格范围覆盖了青藏高原东南部探空站点集中地区。驱动模式的初始条件和边界条件为ERA-Interim再分析资料。模拟时段为2017年7月21日00:00 至25日18:00，其中22—23日青藏高原东南部普遍为晴空，24—25日西藏中南部(即D03的东部)有中等以上强度降水，其余地区(包含D03西部)普遍有小雨，同时考察了晴空和降水条件下的边界层发展。主要物理参数化方案分别为：Dudhia短波辐射方案[53]，RRTM长波辐射方案[54]，Noah陆面过程方案[55]。模拟试验中对比12组不同的微物理、行星边界层和积云参数化方案的组合，确定研究时段内本地适用性最高的参数化方案组合。

图4 模拟区域设置(阴影为地形高度，单位：m)

表3 模式嵌套区域设置

3.2 参数化方案组合对模拟结果的影响

首先检验各方案组合得到的青藏高原边界层内气象要素场的模拟结果(表4)，利用2组探空数据平均值检验的要素为地面、500 hPa、400 hPa温度、露点温度、风向和风速。

表4 模拟时段内边界层12个组合方案中主要8个组合方案的模拟要素在5站点的平均值检验

12种方案组合对边界层内温度模拟的Bias和RMSE均较小，且模拟值偏差随高度升高而减小；温度模拟误差的离散程度对边界层方案的敏感性较显著，其中MYJ

3.3 参数化方案模拟差异的诊断

为了解释模拟结果能较好地反映高原边界层物理过程的机制，对12组参数化方案模拟的近地层感热通量、地气温差和边界层高度进行诊断(图5)。

图5 不同参数化方案组合模拟的D03区域2017年7月地面感热通量(a)、地气温差(b)及边界层高度(c)平均值的时间序列

可以看到，12组参数化方案模拟的各要素变化规律较为一致，但在峰值处有一定差异，这种差异对于认识青藏高原边界层内物理过程并据此修正再分析资料具有重要意义。各方案组合的感热通量模拟结果差异主要体现在06:00的日峰值，其中MYJ边界层方案模拟的感热通量最高；BMJ的感热通量高于浅对流KF结果；两微物理方案的结果相差不大，但是WSM6方案的感热通量下降较Lin方案更加明显。WRF模式中地表感热通量的模拟主要依赖于陆面方案，由于12组试验一致使用Noah方案，因此模拟结果差异不大。由陆面方案向边界层方案传递的地表通量的模拟差异会进一步引起边界层高度的模拟差异。

地气温差对边界层方案的选择较为敏感，YSU模拟值最高，ACM2方案略次之，MYJ明显偏低；而不同云微物理方案和积云对流参数化方案造成的差异很小。边界层方案中局地方案的地表冷却强度要略高于非局地方案，即MYJ方案模拟得到的地表温度相对于YSU方案更低，而夏天2 m气温比地表温度更低，因此，MYJ方案的地气温差比YSU方案显著偏小。对比图2(b)中ERA-Interim计算得到的夏季边界层高度日变化，模拟结果更符合昼高夜低的变化规律，即模拟物理框架普遍可以修正再分析资料描述边界层高度演变的偏差。

对于边界层高度模拟结果，晴空时边界层方案和积云对流方案的作用大于微物理方案，但是微物理方案模拟的边界层高度对降水的响应更加显著。边界层方案中ACM2模拟边界层高度最高，YSU次之，MYJ方案较低；积云对流方案中BMJ模拟结果高于浅对流KF方案；微物理方案中WSM6的响应强于Lin方案。非局地边界层方案YSU和ACM2给出的边界层内湍流混合强度和卷夹强度较强，其模拟的湍流混合强度较MYJ方案更强，因此即使非局地方案中输入的地表感热通量较低，其模拟的湍流混合层仍然相对于MYJ方案更高。WSM6方案中的基本物理过程与Lin方案基本一致，但WSM6方案改进了加热项的计算精度和稳定度，且对于一些经验参数的选择有所不同。此外，虽然两微物理方案的模拟效果均较可靠，但Lin方案对于模拟分辨率的敏感度较高，随着模拟分辨率的提高，WSM6的优势更加明显[56]。

夏季青藏高原对流边界层内湍流发展旺盛，湍流运动发生的动力条件之一是气层间的风速切变，因此，能够较为准确地描述风速切变强度是利用再分析资料分析高原对流边界层的重要基础。由表2可知，模式物理过程对青藏高原东南部边界层内物理量垂直分布的修正作用主要体现在不同层次的水平风场上，考虑物理过程的模式模拟改善了再分析资料中风场部分。对比晴空(7月22日)和降水(7月24日)情形下边界层内3个层次上的水平风场分布(图6)，再分析资料和模拟结果可以刻画由地面到高空风速随高度递增的特征，但风速和风向尚存在较大误差。ERA-Interim再分析资料中青藏高原东南部边界层内各层水平风场风向较为一致，且各层次之间风切变不明显。晴空情形下，400 hPa为一致偏南风，500 hPa白天到夜间风向由偏南风转为偏东风，地面风向日变化小。有降水发生时，400 hPa偏南风夜间强于白天，500 hPa高原中部和东部风向几乎反向。而模拟结果显示边界层内各层次风向不统一，且风速较再分析资料更大，即存在较强湍流运动。晴空情形下400 hPa为一致较强偏东风，对比7月22日06:00和18:00地面风场可以看到风向近乎反向，可以表达出地形影响下的风向变化特征。降水情形下，夜间400 hPa呈现一条明显的风向辐合带。此外，模拟结果各层次之间风速切变也更加明显，青藏高原地区地形等因素引起的山谷环流作用不可忽视。因此在地形作用下30°N—34°N范围内31°N南北的地面风向应接近相反，图6中ERA-Interim的水平风场不能表达这种差异，而模拟结果中却可以看到30°N—34°N范围内31°N以南日间主要为偏南风，以北主要为偏北风。模式物理框架中的边界层内湍流强度和风速切变强于再分析资料ERA-Interim，且与探空观测数据间差异更小。

图6 第三层嵌套区域边界层内不同层次的水平风场分布(左边两列为2017年7月22日代表晴空日，右边两列为2017年7月25日代表降水日)

4 再分析资料修正

模式物理框架对高原东南部边界层内湍流运动的描述相对于再分析资料更加准确，但更多研究中直接利用ERA-Interim数据讨论青藏高原对流边界层发展及其对下游的影响，为得到更加接近实际的结果，可以利用水平风场的最优模拟结果对再分析资料进行修正。利用模式模拟的关键物理量来修正再分析资料的方案，采用一元线性回归方程订正再分析资料中的边界层内水平风场，利用探空资料再次检验并重新计算边界层高度，与原结果进行对比以检验订正效果。选择7月22—25日ACM2+WSM6+BMJ方案模拟的D03区域内650、600、550、500、450、400 hPa水平风场(因变量，记为y)及ERA-Interim数据中的相应值(自变量，记为x)进行计算(表5)，得到的回归方程均通过α=0.05的显著性检验。

表5 ERA-Interim边界层内水平风场与模拟风场间的一元线性回归方程

利用探空资料对研究时段内的拟合风场进行检验(表6)，可以看到除昌都站以外，其余4站的修正后再分析资料对水平风场的描述效果均有明显提升，昌都站的修正效果差可能与一元线性回归方程的拟合效果有限有关。从区域平均值来看，拟合风场对于边界层内风场的修正有效。

表6 模拟时段内拟合风场的探空检验STDE值及其相对于原始再分析资料的减率

利用订正后的再分析资料计算2017年7月22—25日边界层高度见图7。原再分析资料计算的高原边界层高度夜间在500 m左右，白天高值中心可达2000 m(图略)。利用拟合风场数据计算得到边界层高度在整个研究区域内显著高于前者，虽然不一定更接近于真实值，但可以说明模拟物理框架对于风场的约束作用体现在边界层发展高度上。从空间分布来看，夜间高原南部边界的边界层高度梯度极大值修正后减弱，白天高原东北部的梯度加强。修正后的边界层高度普遍在1500 m以上，高值中心位于高原南部；白天青藏高原上空边界层高度有2个大值区，一个在7月22—24日内稳定于高原东北部，最大值达3000 m以上，另一个大值区在7月24日发生降水后向西偏移，即向降水强度更强的高原中部地区偏移，更加符合实际的物理过程。因此，认为通过对行星边界层方案、微物理方案和积云对流方案组合的合理选取，可以修正和改善再分析资料对于高原夏季降水的模拟效果。何由等[57]通过对青藏高原一次强降水过程进行模拟和对比评分发现，WRF模式基本能够重建强降水过程的中心、强度及降水范围，Eta-Ferrier云微物理方案和Betts-Miller-Janjic积云参数化方案的组合在此次模拟降水过程中最优，且行星边界层方案的合理选择能够明显提高降水模拟的效果。这佐证了利用ACM2+BMJ+WSM6方案组合对再分析资料修正有效的结论。

图7 2017年7月22—25日订正后的再分析资料计算的边界层高度(单位：m)

5 总结与讨论

选取ERA-Interim、JRA-55和MERRA-2 3套再分析资料，讨论其在2017年夏季青藏高原边界层发展个例中的适用性，主要结论如下：

(1)在2017年夏季青藏高原东南部边界层内，3套再分析资料对于气象要素的描述能力气温>露点温度>水平风场。考虑到JRA-55露点温度误差较大且MERRA-2提供的垂直层次有限，研究时段内适用性较好的再分析资料为ERA-Interim。

(2)选择3种边界层方案、2种微物理方案和2种积云对流方案进行模拟试验并进行检验，结果表明12种方案模拟水平风场的误差离散程度的改进效果最明显，温度次之，露点温度模拟不理想。因此对于模拟改进的关键物理量水平风场而言，研究时段内本地适用性最高的参数化方案组合是ACM2+WSM6+BMJ。

(3)考虑了青藏高原物理过程的模拟结果修正再分析资料的方案，并利用原始再分析资料线性拟合最优方案组合模拟水平风场，经由部分探空资料检验，改进了青藏高原上边界层发展高度的计算。

利用物理框架对再分析数据约束，可以提高其对边界层水平风场和风在垂直方向上分布的描述能力。然而边界层发展受地表覆盖等下垫面条件影响外，边界层上部的夹卷效应也有重要影响。后续工作中，需要进一步探寻更好的方法来改善再分析资料的露点温度等物理量场，获取更多观测信息来探讨高原气象要素的分布和变化，检验包括CRA-40等再分析资料的适用性。