1 引言
近地面风是重要的气象驱动要素之一， 精细化风场模拟在数值天气预报、 气象灾害风险评估、 大气环境评价等应用中至关重要（Simpson et al， 2016； Yang et al， 2016； Coccia， 2021； 刘馨尹和张宁， 2021； 李艳等， 2015； Hoolohan et al， 2018）。风场在不同地形下变化快速且复杂， 虽然现在已有的观测资料能提供多尺度、 高精度的信息， 但仍难以描述连续的时空分布场（Fernando et al， 2019）。数值模拟作为现在最有潜力的研究方法之一， 可以提供高精度的时空连续场， 在气象要素模拟、 天气预报、 气候分析等众多领域发挥重要作用（王颖等， 2010； 邢雯雯等， 2021； 粟运等， 2022； 张晨炜等， 2022； 王树舟和于恩涛， 2014）。然而， 如何利用数值模拟的方式得到既有湍流细节又与观测相近的高分辨率近地面风场是一个值得研究的问题。另一方面， 随着计算机技术的突飞猛进和数值预报模式的不断完善， 人们对局地精细的气象信息需求越来越高（金荣花等， 2019）， 主要体现在灾害预警的时效性和重大活动气象保障工作上。精细准确的风场模拟研究不仅可以为实况产品提供数据支撑， 更能为精细化网格预报服务奠定基础（师春香等， 2019）。
大气边界层是地面与大气之间物质进行动量、 热量和水汽交换的重要场所（Stull， 1988）。边界层的稳定性和斜压不稳定性决定着风场的变化， 因此在模式中边界层方案对边界层内部状态和流场的模拟精度起着决定性的作用（Carvalho et al， 2012）。中尺度模式中一个网格单元往往代表标量的平均态， 同时也因分辨率太低而无法对三维湍流解析（Mirocha et al， 2010）。大涡模拟（Large Eddies Simulation， LES）是基于大尺度涡旋包含大部分湍流动能， 而小尺度涡旋主要起到耗散作用的原理， 利用空间滤波函数把三维大气中的大尺度涡旋和小尺度涡旋分开， 大涡用湍流运动方程（纳维-斯托克斯方程）显式求解， 小涡则在湍流求解方程中添加应力项求解的一种湍流数值模拟方法。
中尺度模式WRF（Weather Research and Forecast Model）， 一方面可以利用全球预报或再分析资料作为初边值条件， 拥有大气背景场； 另一方面又具备LES所需的动力框架及物理方案， 利于中微尺度耦合（Mirocha et al， 2010）。因此， 前人基于WRF-LES对真实大气环境下已展开了一些探索。 Liu et al（2011）利用WRF-LES对科罗拉多某风电场进行数值模拟， 比较了湍流扩散参数及嵌套方式对风场模拟的影响， 发现双向嵌套得到的湍流细节利于反馈。孙学金等 （2017）利用WRF-LES进行风场模拟， 发现改变土壤湿度会影响近地面风场的模拟精度。刘郁珏等 （2018）针对地形和垂直分辨率进行边界层风场敏感性试验， 结果表明复杂地区需适配同分辨率地形和适当加密垂直分辨率。Hald et al（2019）通过对天气过程进行模拟， 发现WRF-LES在保留中尺度天气背景的条件下， 可以较好地捕捉到小尺度气象特征， 例如山谷环流。
以上研究表明了WRF-LES对局地湍流捕捉以及小尺度天气刻画的优势， 针对局地精细化的风场数值模拟， 必须根据局地实际情况进行最合适的分辨率和参数设置。LES作为WRF中边界层方案的一种， 在1 km这个介于中尺度和微尺度之间的过渡分辨率， 开启大涡模拟方案对里层嵌套模拟效果是否更好还有待进一步研究。同时， 大涡模拟中是包含亚网格通量的， 必须通过亚网格模型对其进行封闭建模， 而不同亚网格方案对近地面风场的模拟效果还有待检验。
鉴于以上， 本文基于WRF-LES进行五层双向嵌套， 选择北京门头沟局地区域的晴天大风个例针对1 km分辨率是否开启LES和千米尺度以下WRF-LES中不同亚网格方案展开近地面风场的模拟， 并利用地面观测进行详细评估检验。本研究主要聚焦两个问题： 其一， 1 km分辨率开启大涡模拟是否对模拟结果有优化； 其二， 大涡模拟中亚网格方案对于近地面风场模拟效果如何， 何种方案对风场模拟结果最佳， 从而加深对WRF-LES的理解并为其应用提供科学的参考依据。
2 模拟区域与观测数据
本文的研究区域位于北京门头沟区东部， 地处华北平原向蒙古高原过渡地带， 地势西北高， 东南低。该地区属中纬度大陆性季风气候， 高空大气环流属于盛行西风带， 冬季干燥多风， 夏季潮湿多雨。门头沟冬春过渡之季阳光明媚， 多伴大风天气。选择门头沟区域的原因如下： 其一， 冬春换季常出现大风天气， 干燥大风天气对于模式模拟具有代表性。其二， 该区域观测密集， 验证具有可靠性。其三， 该地区有平原， 丘陵和山地等多种地形， 利于检验模式在不同地形下对风场的模拟能力。本研究使用中国分钟级地面自动气象站的10 m风速风向数据， 数据的原始采样间隔为1 min， 在结果验证时采样成与其对应的5 min分辨率。根据不同海拔将站点进行分类， 具体观测站点信息如表1 。
表1 站点基本信息Table 1 Basic infomaition of observations |
| 站号 | 经度/°E | 纬度/°N | 海拔/m | 分类 |
|---|---|---|---|---|
| A0358 | 116°6′14.04″ | 39°56′27.96″ | 98 | 平原 |
| A1381 | 116°6′33.84″ | 39°54′51.12″ | 96 | 平原 |
| 54505 | 116°9′23.04″ | 39°53′16.08″ | 86 | 平原 |
| A1036 | 116°9′27.00″ | 39°55′59.88″ | 91 | 平原 |
| A1066 | 116°6′50.04″ | 39°51′25.92″ | 119 | 山麓 |
| A1109 | 116°7′32.16″ | 39°57′47.16″ | 143 | 山麓 |
| A1067 | 116°4′18.12″ | 39°52′08.04″ | 427 | 丘陵 |
A为区域站， 54505为国家一般气象站 |
3 模式设置和试验设计
3.1 模式参数设置
本研究采用WRF-LES（WRF版本3.9.1.1）对门头沟区域进行五层嵌套， 通过动力降尺度至最里层为111 m， 其中中尺度区域（千米尺度）为D01， D02， D03三层嵌套， 微尺度（百米尺度）为D04， D05两层嵌套。初始场和侧边界条件来自国家气象信息中心研制的全球大气实况分析产品， 该产品是基于四维集合变分混合同化方法， 同化实时收集的中国和全球常规观测资料和卫星探测数据后， 生成的覆盖全球的三维大气实况产品（Liang et al， 2020）。该实况产品的空间分辨率为13 km×13 km， 时间分辨率为6 h。
模式模拟区域如图1 所示， 模式中心位置为39°54′ 51.12″N， 116°6′33.84″ E， 其中第一层包含整个河北省， 第二层包含整个北京市， 最里层位于门头沟区新城区， 北邻龙泉镇， 南接潭柘寺镇。Wyngaard （2004）认为各嵌套层水平网格数一致利于大尺度湍流的形成和发展， 因此本研究网格的五重嵌套水平网格数量均为121×121， 垂直设置40层。每一层设置的物理方案如下表2 ， 在中尺度嵌套层（D01~D03）使用的物理参数方案如下表格， 在千米尺度以下的嵌套层（D04~D05）关闭边界层， 开启大涡模拟， 时间集成使用3阶龙格-库塔方案（Wicker and Skamarock， 2002）， 离散方式使用显式六阶离散（Knievel et al， 2007）， 使用到的亚网格方案具体见3.2。Liu et al （2011）发现双向反馈有利于大涡模拟对于涡旋细节的捕捉， 故本研究选择使用双向反馈。
图1 模式嵌套示意图（a）和最里层D05的地形分布（b）空心三角形代表区域站， 实心三角形代表国家一般气象站 Fig.1 Simulation domains （a） and the terrain distribution of D05 （b）.Hollow triangle represents the Regional Station， solid triangle represents the National General Weather Station |
表2 模式逐层嵌套网格参数设置Table 2 Domain size， grid spacing and time steps used in the WRF model simulations of each nested domain |
| 嵌套层 | 网格数 | 水平分辨率 | 水平跨越/km | 时间步长/s |
|---|---|---|---|---|
| D01 | 121×121×40 | 9 km | 1089 | 18 |
| D02 | 121×121×40 | 3 km | 363 | 6 |
| D03 | 121×121×40 | 1 km | 121 | 2 |
| D04 | 121×121×40 | 333 m | 40.29 | 0.667 |
| D05 | 121×121×40 | 111 m | 13.43 | 0.222 |
表3 WRF和WRF-LES模式物理参数设置Table 3 Physical schemes used in WRF and WRF-LES models |
| 模式 | 物理参数化方案 |
|---|---|
| WRF | Kain-Fritsch积云方案 |
| YSU边界层方案 | |
| Thompson微物理方案 | |
| RRTM长、 短波辐射方案 | |
| MM5近地层方案 | |
| Noah陆面过程方案 | |
| WRF-LES | 3阶Runge-Kutta时间集成 |
| 1.5TKE， SMAG， 1.5TKE_NBA， SMAG_NBA湍流方案 | |
| 6阶数值离散方案 |
3.2 亚网格模型介绍
3.2.1 标准的WRF-LES亚网格模型
WRF自带的基本亚网格（Subgrid-Scale， SGS）应力模型有两种， 其一是Smagorinsky模型（Lilly， 1967； Smagorinsky， 1963）， 其二是基于1.5阶TKE的模型（Lilly， 1967）。两者基于亚网格应力的涡旋方程如下：
式中： 代表亚网格应力； 是涡黏系数； 是解析的应变速率张量。在Smagorinsky模型中， 涡黏系数使用以下方程：
式中： 是Smagorinsky系数， 在WRF中这个参数默认0.25； 是湍流长度尺度， ； Pr是数值为0.7的湍流普朗特数； 是Brunt-Vaïs la频率。
1.5阶TKE的涡黏系数方程如下：
式中： e是亚网格的湍流动能， =0.15， 当选择这个子过滤器的选项时， WRF会集成一个亚网格湍流预测方程。
两个模型都是易于理解和实现的， 其中Smagorinsky模型相对于1.5阶TKE模型应用更加广泛。但是两者仍有局限， 主要体现在以下这几个方面： 1）它们使用的是常量， 在普遍情况下是最优； 2）两种模型都假定湍流耗散仅在局部平衡下产生； 3）没有包括后向散射（Mirocha et al， 2010）。
3.2.2 NBA模型
基于第三个缺点， 考虑后向散射的亚网格模型嵌入至WRF中， 该模型实现了非线性后向散射和各向异性（Nonlinear Backscatter and Anisotropy， NBA）的应力模型（Kosović， 1997）， 该模型基于非线性本构关系， 并包括考虑剪切均匀湍流情况下的后向散射和法向应力的二阶项。NBA模型用亚网格湍流动能表示为：
式中： 是可解的旋转张量； ； ； ； 。该模型实际上只依赖于后向散射系数 ， 在该系数无限逼近于0时， 由方程式给出的NBA模型（4）和（5）分别回归到线性涡粘Smagorinsky模型和1.5阶TKE亚网格应力模型。后向散射参数 代表的是平均值， 它降低了整体耗散率， 导致近表面应力和整体光谱特征的改善（Kosović and Curry， 2000）。在对流、 中性和稳定分层的理想大气边界层模拟中， NBA模型比线性涡度模型有显著改进， 包括与近地表的相似性剖面有更好的一致性， 改进了频谱特征和各向异性（Mirocha et al， 2010）。
3.3 试验设计
在WRF的边界层方案中， YSU是一阶闭合的非局地方案， 利用非局地的涡旋扩散方法去显式求解边界层的顶部， Hu et al （2013）认为， YSU通过垂直混合从而改善了稳定状态下的风速估计， 是专注于改善风力的PBL方案之一。众多学者已经证明YSU边界层方案的优势（Carvalho et al， 2014a； Gholami et al， 2021； 张少济等， 2017； 李玉鹏等， 2018； 周荣卫等， 2018）。1 km是介于中尺度和微尺度之间的过渡范围， 作为两个尺度之间的桥梁， 它的精准模拟决定了最里层嵌套的模拟效果。同时， 在1 km这个尺度是有多重选择的， 既可以使用普通边界层方案， 也可以选择大涡模拟［使用手册中建议在2 km以下分辨率使用（Skamarock et al， 2008）］。大涡模拟的处理并非一个必要的物理选择， 所以它在中小尺度过渡模拟的适用性有必要研究。因此， 在1 km分辨率的嵌套层进行大涡模拟的敏感性试验， 为对比1 km打开动力模块与否对风场模拟的影响， 设计试验1 km-YSU和1 km-LES。前者在第三层嵌套使用YSU边界层方案， 关闭大涡模拟； 后者在第三层嵌套关闭YSU边界层方案， 开启大涡模拟。嵌套的最里层（D04， D05）默认都关闭边界层， 开启动力模块。
大涡模拟的优点在于筛选出对动量、 质量和能量影响大的大尺度涡旋， 过滤掉主要起耗散作用的小涡。过滤的小涡通过大涡流场运动方程中附加的应力项体现， 这个应力项就是通过亚网格模型来建立。因此， 使用何种亚网格模型解析小涡对大涡的运动方程求解也至关重要。本研究选择的四种SGS模型有： 1.5阶湍流闭合模型（1.5TKE）， Smagorinsky一阶闭合模型（SMAG）和非线性后向散射各异向模型的两种形式（1.5TKE_NBA， SMAG_NBA）。
表4 边界层敏感性试验设计Table 4 Experimental design of atmospheric boundary layer sensitivity |
| 组别 | 试验名称 | 边界层参数 | ||
|---|---|---|---|---|
| 中尺度嵌套 （dx， dy=9 km， 3 km） | 1 km嵌套 （dx， dy=1 km） | 微尺度嵌套 （dx， dy=333 m， 111 m） | ||
| 试验一 | 1 km-YSU | YSU | YSU | LES（1.5TKE） |
| 1 km-LES | YSU | LES（1.5TKE） | LES（1.5TKE） | |
| 试验二 | 1.5TKE | YSU | YSU | LES（1.5TKE） |
| SMAG | YSU | YSU | LES（SMAG） | |
| 1.5TKE_NBA | YSU | YSU | LES（1.5TKE_NBA） | |
| SMAG_NBA | YSU | YSU | LES（SMAG_NBA） | |
4 模拟结果检验与分析
4.1 1 km边界层方案敏感性试验
1 km作为中尺度和微尺度的桥梁， 对其选择合适的参数化方案有利于对微尺度模拟的反馈。第一个试验是针对模拟的第三层（1 km）是否打开大涡模拟而进行的， 该试验主要是为了评估WRF-LES在1 km尺度的适用性。首先， 为评估外层嵌套区域作为里层的初始边界条件的准确性， 利用D01~D03区域逐小时的风向风速进行检验。图2 是根据所有逐小时的模拟结果（D01~D03）与站点的逐小时的风速进行统计得出， 可以看到， 第三层是否打开大涡模拟， 会对第一层和第二层嵌套产生影响， 这是因为WRF模拟中使用的是双向反馈， 其中第三层使用YSU边界层方案整体的模拟结果与观测更为接近。相比于开启大涡模拟， 第三层关闭大涡模拟的归一化均方根误差更低， 相关系数也更高。对于归一化标准差而言， 第三层开启大涡模拟的标准差较高， 说明整体的模拟风速波动性大， 这符合大涡模拟对风速脉动的刻画。
为检验两者描述逐小时风向的变化能力， 分别从平原， 山麓和丘陵选择具有代表性的站点： 54505、 A1066和A1067。图3 为不同地形站点的风矢量时间序列图。从图3 中可以看出， 针对1 km是否打开大涡， 两者在D01~D03的风向分布大体一致。在平坦地区［图3 （a）， （b）， （d）， （e）］， 1 km-YSU的D03在风速和风向上都更接近观测， 1 km-LES则在风速和风向上偏离程度较大。在较复杂的地区［图5 （c）， （f）］， 1 km-LES在风向上的波动性明显增加， 在16日20:00（北京时， 下同）至17日04:00的风速风向上变化捕捉能力比1 km-YSU更佳。在17日00:00 - 08:00， 风速较小的时刻， 1 km-YSU出现风速高估现象。由此可见， 1 km-YSU在平坦地形表现较好， 1 km-LES在复杂地形表现更佳。
图3 站点54505（a， d）、 站点A1066（b， e）和站点A1067（c， f）的逐小时风矢量随时间分布橙色代表观测， 深蓝色代表不同嵌套层的模拟结果 Fig.3 Time series of wind directions from simulation and observations at 54505（a， d）， A1066（b， e） and A1067（c， f）.Orange represents observed data， darkblue represents simulation results from different domain |
图4 1 km-YSU（a~b）中D04和D05与1 km-LES（c~d）中D04和D05的风矢量散点分布（单位： m·s-1）Fig.4 1 km-YSU（a~b）D04 and D05 and 1 km-LES （c~d） D04 and D05 Wind vector scatter diagram.Unit： m·s-1 |
图5 1 km-YSU（a）和1 km-LES（b）的近地面风速偏差和均方根误差空间分布（单位： m·s-1）红色为正偏差， 蓝色为负偏差， 圆圈中数值代表均方根误差 Fig.5 Spatial distribution of near-surface wind deviation and root mean square error of D05 in （a）1 km-YSU （b）1 km-LES.Unit： m·s-1.Red circle means the positive deviation， blue circle means the negative deviation， and the value in the circle represents the RMSE |
其次， 为评估最里层百米风场的模拟结果， 分别从整体评估指标， UV矢量分布， 误差空间分布和平均风速日变化去检验。表5 给出1 km边界层敏感性试验整体评估结果， 平均偏差（MB）、 均方根误差（RMSE）和平均绝对偏差（MAE）三者计算公式如下：
式中： 为模式模拟值； 为观测值； 为样本数。从表5 中可以看出， U分量呈现高估， V分量呈现低估， 1 km-YSU在平均绝对偏差和均方根误差与1 km-LES相差不大， 但数值上较1 km-LES要低。图4 为百米尺度的嵌套层的UV矢量分布。由图4 可知， 1 km-YSU和1 km-LES在D04和D05的分布变化差别不大， 整体为偏北风， 以西北风为主。大体上看， 两者的矢量分布主要集中在西北方向， 西南和东南方向。其中1 km-LES在西南方向矢量分布较1 km-YSU更为分散， 与观测更为吻合。但是两者在大风（>15 m·s-1）以及东北方向矢量的捕捉能力都不佳。在千米尺度以下， D04（333 m）和D05（111 m）模拟结果差别不大。两者相似的原因可以归结于， 双向嵌套形式使得内外层可以传递信息并相互反馈， 父域和子域之间会有更高的一致性， 从而减少了从域边界反射产生的错误（Harris and Durran， 2010）。综上， 1 km-YSU风矢量在主导风（西北风）和次主导风（东南风）更加集中， 而1 km-LES相对而言较为分散。
表5 1 km边界层方案试验D05近地面风速统计检验结果Table 5 Statistic results of different 1 km boundary layer schemes simulated near-surface wind in D05 |
| 风速 /（m·s-1） | 试验名称 | MB | MAE | RMSE |
|---|---|---|---|---|
| U分量 | 1 km-YSU | 0.8312 | 3.0109 | 3.8487 |
| 1 km-LES | 1.2825 | 3.3975 | 4.2218 | |
| V分量 | 1 km-YSU | -1.6566 | 3.3129 | 4.2202 |
| 1 km-LES | -1.49 | 3.3604 | 4.2398 | |
| 全风速 | 1 km-YSU | 1.2145 | 2.8428 | 3.616 |
| 1 km-LES | 1.3801 | 2.9388 | 3.7067 |
粗体表示多组实验中最优试验 |
不同方案的偏差和均方根误差空间分布如图5 所示， 红色为正偏差， 蓝色为负偏差， 圆圈数值代表均方根误差， 可见空间上两者分布大体相同。在两者误差分布上， A1109， A0358站点都出现低估， A1103， A1022站点都是高估最为明显。针对站点的评估， 两者差异较小， 但1 km-YSU表现更优， 体现在偏差和均方根误差数值更小。总体而言， 1 km-LES的均方根误差要高于1 km-YSU， 原因可能是关闭边界参数， 使得LES无法较好的再现当时天气状况。
通过D05的5 min输出平均风速的日变化（图6 ）可以看出， 模式大体上能够较好的刻画湍流的脉动特征。1 km-LES相较于1 km-YSU高估现象出现频次更多， 在20:00至04:00两者高估现象明显， 但1 km-YSU的风速波动情况更靠近观测。模式对风速下降的趋势是能够捕捉到的， 但在时间上出现了延后一小时的现象。总体而言， 模式对于5 min风速脉动捕捉较好， 在保留一定脉动细节的基础上， 1 km-YSU比1 km-LES更接近观测分布。
综上所述， 从整体统计， 空间分布和时间序列都可以看出1 km-YSU比1 km-LES模拟效果要更好， 尤其是在平坦地形。
4.2 亚网格模型敏感性试验
不同的亚网格方案会因为解析湍流方式的不同从而导致对风场的影响。不同亚网格方案在最里层嵌套的逐5 min风速泰勒图（图7 ）分布可见四个模型对近地面风速模拟的影响差异较小， 1.5TKE和SMAG相近， 两个NBA方案相近。就归一化标准差而言， 两个NBA方案较高， 这是因为NBA对后向散射能量的解析， 使其捕捉到更多的脉动细节， 风速波动大而使得误差增加。从归一化均方根误差可以看出， 1.5TKE和SMAG分布相近， NBA两个模型较前两者误差略高。从相关系数来看， 三者差别不大。
表6 不同亚网格方案D05近地面风速统计检验结果Table 6 Statistic results of different SGS models simulated near-surface wind of D05 |
| 风速 /（m·s-1） | 试验名称 | MB | MAE | RMSE |
|---|---|---|---|---|
| U分量 | 1.5TKE | -1.6566 | 3.3129 | 4.2202 |
| SMAG | -1.1333 | 3.1593 | 4.0579 | |
| 1.5TKE_NBA | -2.0655 | 3.4245 | 4.2763 | |
| SMAG_NBA | -2.0655 | 3.4148 | 4.2451 | |
| V分量 | 1.5TKE | 0.8312 | 3.0109 | 3.8487 |
| SMAG | 0.5117 | 2.9424 | 3.8355 | |
| 1.5TKE_NBA | 1.1544 | 3.0056 | 3.7276 | |
| SMAG_NBA | 1.1486 | 3.0996 | 3.8361 | |
| 全风速 | 1.5TKE | 1.2145 | 2.8428 | 3.616 |
| SMAG | 0.6936 | 2.639 | 3.3505 | |
| 1.5TKE_NBA | 1.518 | 3.0763 | 3.8515 | |
| SMAG_NBA | 1.5123 | 3.1099 | 3.8626 |
粗体表示多组实验中最优试验 |
风矢量分布（图8 ）中， 条带距离原点的距离为该方向（某一个中心角度左右22.5°， 即夹角为45°）上的平均风速。带宽表示该角度风向出现的频率， 颜色的深浅代表该风向上出现风速大小的一致性［参考meteva说明文档（https：//www.showdoc.com.cn/meteva）： 某风向的范围内的所有样本的风速的平均值m和标准差d， 采用 m/（m+d）来代表风速的一致性］。结合表6 与不同的亚网格方案的风矢量分布统计图可以看出， 四者都能较好地模拟出主导方向（西北）的大风， 西北方向平均风速大小， 出现频次与观测最为接近的是1.5TKE， SMAG略有低估， 而两个NBA略有高估。对于次主导风（东北风）， 四者均出现高估现象， 1.5TKE_NBA的几乎与观测一致。针对偏南风（南风， 东南风， 西南风）， 四者都出现低估， SMAG模拟效果最佳， 但1.5TKE_NBA的一致性最高。总体而言， SMAG表现最佳， 出现的高估和低估现象都最小， 与观测最为吻合。
风矢量特征分析图只能看出风速风向模拟的整体情况， 不能显示风速随着时间变化。从图9 可以看到， 风速误差随时间的变化， 其中误差棒圆点为每三小时内所有站点误差的中位数， 上下两头分别为75和25分位数。结合误差时间序列图（图9 ）和平均风速时间序列图（见图6 ）可以看到， 风速较大的时刻（14:00 -16:00）， 模式误差分布范围大， 误差中位数低于0， 大体表现为低估； 风速较小的时刻（23:00至次日04:00）， 误差分布范围较小， 同时25分位数都高于0， 说明高估现象明显， 这一点与前人的模拟结果相似（Doubrawa et al， 2018）。SMAG在所有的亚网格模型中表现最佳， 误差分布范围相较而言最小， 中位数分布大部分都接近于0。
图9 不同方案的误差棒日变化误差棒圆点为每3 h内所有站点误差的中位数， 上下两头分别为75和25分位数 Fig.9 Diurnal variation of 10-m wind speed errors （m·s-1） of each OBS with different SGS models at D05.Markers represent the median error， and whiskers give 25th and75th percentile values for each bin.Bars indicate how much of the total data considered belong in each bin |
由风速风向的概率密度分布（图10 ）可知， 风速主要集中在0~5 m·s-1之间， SMAG与观测最为接近。模式在9 m·s-1和15 m·s-1左右各出现小高峰， 这说明模式对于中高风和强风出现高估现象， 这与前人研究的结果类似（Carvalho et al， 2014b； Hald et al， 2019； Giannaros et al， 2017）。
从不同方向的风速瞬时分布来看， 不同的亚网格模型对于风场的刻画有一定差异。由2月17日04:20（夜间的稳定大气条件下）的风场经向、 纬向和垂直分量分布图（图11 ）可知， 1.5TKE和SMAG与其对应的NBA模型各方向风场趋势大体一致， 但在细节上略有不同。就U风量而言， 都呈现出山区西风分量， 平原东风分量的特征， 西部区域1.5TKE和SMAG模拟中是大块的西风分量（正U）， 而NBA中明显看到西风分量被打碎。西北侧两座山头之间形成的狭管地形四种方案都有一股东风分量， 顺着山坡滑泄在局地形成环流。SMAG_NBA的东风分量分布呈现为三角形， 与扇形的山谷地形分布更加符合。就V风分量而言， 在东南侧平原地区， 两个NBA模型在主导的南风分量中， 有细小的北风分量， 平原中的湍流已经受地形影响较小了， 因此几乎是SGS模型对湍流解析的影响， 可见后向散射模型的优势。垂直方向上， SMAG_NBA相比于SMAG湍流解析能力更强， 在保持原有细节的基础上， 把湍流分得更加细碎， SMAG_NBA在水平和垂直方向上对于山谷空气流动的捕捉更为精细。由于本研究重点在于研究近地面风场， 垂直风场仅作辅助信息， 可见SMAG_NBA模型在垂直方向上显现的优势， 真实性还需观测进一步验证。综上， SMAG在数值检验和风场分布来看， 比1.5TKE更接近观测， 对风场的刻画更为精细。因此， SMAG是局地风场模拟的最好选择。
5 结论与讨论
本文利用中尺度WRF内层嵌套微尺度大涡模拟， 通过五层嵌套将分辨率从9 km降到了111 m， 对门头沟东部区域进行了近地面百米风场的模拟。针对1 km分辨率开启大涡模拟是否对风场模拟结果有优化和WRF-LES中亚网格方案对风场模拟的影响两个角度展开试验， 通过近地面观测对模拟结果进行检验， 得出主要结论如下：
（1） 在中微尺度过渡分辨率1 km开启YSU边界层较LES边界方案模拟效果更佳， 在保持捕捉湍流细节的能力下， 拥有较低的误差， 建议在1 km模拟时使用YSU。原因有两点： 其一， YSU边界层可以对地形进行校正， 双向嵌套有利于反馈到里层。其二， LES需关闭边界层， 失去边界物理方案的约束， 解析湍流而带来的误差较大， 但不可否认其对湍流的捕捉能力。
（2） 从整体检验的统计特征来看， 几种亚网格模型差异较小， SMAG亚网格模型表现最佳， 体现为均方根误差最小， 相关系数最高。SMAG的涡黏模型中应力大小是随着静态稳定性而适当进行调整的， 湍流会有一定程度的限制， 使得误差和波动都会相对于1.5TKE要小。鉴于此， SMAG模拟时所耗费的时间会相对于1.5TKE要长一些， 在模拟不稳定的对流情况需根据实际而选择。
（3） WRF-LES依然存在对风速模拟的常见现象， 即高值低估， 低值高估， 四种亚网格方案都存在这个问题。其中， SMAG相较之下与观测更接近， 在风矢量分布、 风速风向概率分布和风场瞬时空间分布都表现较好， 而 NBA模型高估现象较为突出。但两个NBA方案都表现出对细节捕捉的优势， 在主导风向的风速模拟也与观测最为接近。其精致的细节刻画能力主要来自于NBA模型中新增的后向散射参数。
综上结论， WRF-LES中SMAG相对于对于微尺度真实大气下风场模拟具有较好的优越性， 可以为局地区域进行环流分析， 模拟预报等提供可靠的数据支持。未来可以面向局地风电场， 进行不同高度的风场分布， 湍流变化的模拟评估检验， WRF-LES进行风能资源评估的应用有待开发。现已看出大涡模拟对近地面风场的模拟的优势， 未来可以从风电场的角度去（风切变/湍流动能变化/垂直风速廓线）检验分析WRF-LES的适用性。另外， 仅依靠地面观测检验大涡模拟的湍流细节的捕捉较为单一， 未来可以结合测风激光雷达和涡动相关系统利用垂直速度、 热通量和湍流动能等物理量进行模拟检验， 通过湍流特性来详细分析亚网格方案造成差异的机理。