边界层大气受地球表面的影响，对表面强迫力的响应时间约为1 小时或更短，对人类生产生活有直接影响^[1].针对大气边界层的研究不仅有助于加深对边界层物理过程的认识，而且是城市空气质量研究的基础.近年来城市化的推进不仅带来城市热岛等气象问题，也加剧大气污染等环境问题.大气运动直接决定或影响污染物扩散输送过程，如水平、垂直输送，湍流混合，干湿沉降等.在数值模拟中，输入空气质量模式的气象场数据包含显著的不确定性，影响最终的模拟结果^[2].重度污染过程多发生在相对弱的天气形势下(小风，无降水，低边界层高度)，多受大尺度高压控制，气压梯度力小^[3].在这种天气背景下湍流运动和平均运动处于一个量级，因此风向容易发生变化，增加数值模拟的难度.同时高污染浓度下风场较小的模拟误差也会给最终的污染物浓度模拟结果带来较大误差.所以模拟弱天气形势下的重度污染过程是一个气象上的难点.

边界层湍流有利于污染物在边界层内扩散，因此边界层高度决定了污染物扩散的有效空气体积. 受边界层结构日变化影响，城市大气污染浓度呈现明显的日变化特征.夜间边界层高度较低，逆温层发生频率较高，不利于污染物扩散，往往出现重污染过程^[4].姜大膀等人研究发现兰州市区低空逆温强度和污染物浓度间存在显著正相关^[5].白天随着地表辐射加热，边界层高度逐渐抬升，逆温层强度减弱甚至转变为不稳定层结，有利于污染物扩散，从而降低浓度.朱燕舞等^[6]分析北京城市边界层内大气污染物的垂直廓线，发现一天中不同时刻对应不同的污染物垂直分布特征.冬季气温低，太阳辐射弱，污染物的化学转换速率减慢，边界层内的物理扩散输送过程对污染物浓度变化起更大作用^[7-9].因此，准确地模拟边界层结构对大气污染过程的模拟至关重要.

随着数值模式发展和计算机运算能力提升，数值模式在边界层结构研究中得到广泛运用^[10-12].由于气象模式分辨率不足以解析湍流运动，需要引入边界层参数化方案或称为湍流参数化方案，用于计算边界层内热量、动量等物理量的湍流脉动通量^[13].不同边界层方案模拟的气象场差别较大. Zhang等^[14]考察MM5模拟的地表风温场日变化对5种边界层方案的敏感性，发现在对流边界层中，模拟结果很大程度依赖于边界层方案的选择.Hu 等^[15]分析WRF中3种边界层方案对模拟结果的影响发现，由于不同方案模拟的边界层顶卷挟混合强度不同，导致边界层内温度、湿度模拟结果的差异.

沈阳位于我国东北平原，拥有720 万人口，是我国传统的重工业基地，发展金属冶炼、化工、建材等高耗能产业，每年消耗800 万吨燃煤，其中20% 用于采暖和生活^[16]，环境空气质量较差.尤其在冬天采暖季，由于燃煤排放大量SO₂，同时低层大气出现逆温层结的频率较高，较其他季节更容易形成严重污染^[17].因此有必要针对这一区域的冬季边界层特征展开研究.

本文讨论新一代中尺度气象模式WRF 中两种不同边界层方案(局地MYJ方案和非局地YSU 方案)对沈阳地区冬季边界层模拟的影响，主要分析对污染物扩散起关键作用的温度层结、低层风场、边界层高度等要素，试图为空气质量模式提供更加准确的气象场.

边界层方案可分为局地方案和非局地方案:前者假设每个格点上的脉动通量完全由该格点上物理量的平均量决定；后者综合考虑该格点及周边格点对脉动通量的影响.对流边界层中大尺度湍涡起主要作用，局地方案的假设有较大偏差^[18].本文考察WRF模式中两种边界层方案对模拟结果的影响: 局地MYJ方案和非局地YSU 方案.这两种方案均是WRF中较早采用的方案^[19]，在前人研究中得到广泛应用^[20-23]，其中YSU 方案进行了一个重要改进^[24].

MYJ方案采用Mellor和Yamada^[25]的湍流闭合方法表示表面层以上的湍流，湍流扩散系数通过湍流动能计算，边界层高度由湍流动能廓线决定.该方案适用于所有稳定条件和弱不稳定条件的边界层，但在对流边界层中误差较大.

YSU 方案是非局地K 理论方案，通过在控制方程中加入逆梯度项表示非局地通量^[24].YSU 方案是对MRF方案的改进^[26].相比MRF 方案，YSU 方案中热对流产生的混合层高度升高，风剪切产生的混合层高度降低；逆梯度项值的减小，使得边界层结构更接近中性，解决了MRF 方案中由于逆梯度项过大导致层结过于稳定的问题.

WRF模式由1°×1°经纬度分辨率的NCEP/ NCAR再分析资料提供初始场和边条件，设置3重双向嵌套的模拟区域，水平分辨率分别为50，10，2km. 如图 1所示，最外层区域覆盖中国东北，朝鲜半岛，渤海，黄海等区域；最内层区域东西140km,南北200km,包括沈阳及周边铁岭、抚顺、本溪、辽阳、鞍山等城市.垂直方向分27层，模式顶高约20km,采用上疏下密分层以更精确地模拟边界层结构^[14]，其中1000 m 以下分7 层，最低层高50 m.50km 和10km分辨率的区域使用Kain^[27]改进的积云对流参数化方案，最里层区域的水平分辨率能够解析积云对流过程，所以不需要进行参数化.其他物理过程方案设置如下:微物理过程采用Lin 方案^[28]；长波辐射RRTM 方案^[29]；短波辐射Dudhia方案^[30]，每10分钟调用一次；表面层方案需要和边界层方案配合使用，对MYJ方案使用Eta表面层方案^[31]，YSU 方案对应MM5表面层方案^[32]；陆面过程采用Noah 方案^[33].使用美国地质调查局(USGS)提供的24类地表植被和土壤分布等下垫面资料.结合天气条件和污染数据，选取2004年1月2日至16日进行模拟分析.期间无明显降水过程，沈阳地区SO₂ 浓度平均值达到0.13 mg·m^-3，超过国家环境质量标准规定的日均浓度一级标准0.05 mg·m^-3，接近二级标准0.15 mg·m^-3.本文主要关注边界层结构的日变化，相关的分析均基于模拟时段内的平均状态.

图 1

Fig. 1

图 1 WRF模拟区域设置(a)最外层表示分辨率为50km 的区域，里层的两个长方形分别表示分辨率为10km 和2km 的区域； (b)WRF最里层区域，圆点表示12个自动气象站分布，中间的方框表示沈阳城市20km×20km 区域. Fig. 1 Spatial coverage of the domains for the WRF simulation (a) The outer rectangle represents the coarse domain at resolution of 50 km, and the inner rectangles represent the nested domains at resolutions of 10 km and 2 km, respectively. (b) The innermost domain with automatcc weather stations presented by the filled circles. The inner rectangle represents 20 km×20 km Shenyang City area.

为了检验WRF模拟气象场的准确性以及边界层方案的影响，利用沈阳及周边地区12个自动气象站(如图 1b所示)的观测数据对模拟结果进行统计检验，参加检验的气象数据包括2 m 高度温度、 10m高度风速风向.采用相关系数(R)，平均偏差(MB)，平均误差(ME)，命中率(Hit Rate,HR)以及Willmott^[34]提出的相符指数(d)等统计量来衡量WRF的模拟能力.命中率HR 表示在所有参加统计的数据对中，观测和模拟的差值小于某一标准的数据对所占百分比.我们规定标准差值如下:温度2K，风速1m·s^-1，风向30°.相符指数d的范围从0到1，0表示模拟值和观测值完全不相符；1表示模拟值完全符合观测值.各统计量的定义如下所示:

(1)

(2)

(3)

(4)

其中P_i，O_i分别表示模拟值和观测值，P，O分别表示模拟值和观测值的平均值.

统计结果如表 1 所示，WRF 模拟和观测较一致，其中对温度的模拟最准确，但模拟值偏低.观测表明城市化导致城市气温升高^[35-36]，WRF 中使用的土地利用数据反映20世纪80年代的城市分布情况，没有考虑新近城市化的作用，导致模拟温度偏低.另一方面，有研究表明MYJ方案和YSU 方案存在模拟温度偏低的现象^[20-22].Hu 等人^[15]比较WRF中MYJ方案和YSU 方案得到，两者模拟2m 温度的MB分别为-1.25 ℃和-0.9 ℃，均高于本文的结果，即本文的模拟温度偏低程度更大，进一步证实土地利用数据滞后的影响.相比MYJ 方案，YSU 方案温度模拟的MB 和ME 更接近零值，HR 和d更大，模拟更加准确.YSU 方案模拟的温度高于MYJ方案，Hu等人^[15]的研究得到相同结论.这可能是由于YSU 方案中垂直混合更强，位温较高的高层大气和低层大气混合，导致低层气温升高，与Hong和Pan^[26]，Srinivas等人^[37]的模拟结果一致. Kotroni等人^[38]针对MM5中边界层方案的敏感性试验表明非局地方案模拟的2 m 温度高于局地方案且更接近观测值.WRF 风速模拟准确性低于温度模拟，模拟风速偏大，相关系数略低于Kwun 等人^[39]的结果(MYJ:0.6483，YSU:0.6304).与温度模拟类似，WRF 低估了城市对风场的摩擦减弱作用.另一方面，WRF 模式对低层风速的模拟可能存在系统性偏差.Shimada等人^[23]利用WRF 中7 种边界层方案进行为期一个月的风速模拟，发现所有方案模拟的边界层低层风速均偏大.不同边界层方案对风速模拟结果影响不大.风向模拟和观测平均有45°左右的夹角，不同方案的影响不明显，和前人的敏感性试验结论一致^{[40, 19]}.

表 1(Table 1)

表 1 WRF不同边界层方案(MYJ和YSU)模拟气象场和12个自动站观测数据比对统计检验结果 Table 1 Statistic results for temperature at 2 m height (T2), wind speed (WS) and direction (WD) at 10 m height between simulation of WRF with MYJ and YSU boundary layer schemes，respectively，and observation

表 1 WRF不同边界层方案(MYJ和YSU)模拟气象场和12个自动站观测数据比对统计检验结果 Table 1 Statistic results for temperature at 2 m height (T2), wind speed (WS) and direction (WD) at 10 m height between simulation of WRF with MYJ and YSU boundary layer schemes，respectively，and observation

以沈阳自动气象站为例，图 2 给出观测和模拟的气温、风速日变化，可以看到模拟的2 m 气温偏低，YSU模拟值高于MYJ，两种方案间的最大温差接近2 ℃，出现在日出前后的0700LST.中午边界层发展旺盛，向下的动量输送增强，观测表明低层风速增大，模拟也得到相同的趋势，但模拟风速偏大，夜间偏差大于白天，MYJ方案和YSU 方案间的差别不明显.Zhang 和Zheng ^[14]利用MM5 模拟发现模拟风速夜间偏大，白天偏小.这种系统性偏差可能是数值模式固有的问题.

图 2

Fig. 2

图 2 观测和模拟的气温风速日变化 Fig. 2 Diurnal variation of the observed and simulated temperature and wind speed

综上所述，WRF 模拟的气象场可以基本代表真实大气状态，不同边界层方案对温度场模拟影响较大，YSU 方案模拟的低层气温高于MYJ方案，风速风向模拟对边界层方案的敏感性较低.同时注意到由于模式使用的土地利用数据年代较早，未能体现近年来城市化的影响，导致模拟温度偏低，风速偏大.

温度层结决定大气稳定度，影响污染物在边界层内的扩散，低层逆温往往导致重度污染^[4-5].下面分析不同边界层方案对温度层结模拟的影响.

为减小局地扰动，如地形、土地利用的不均匀分布，对模拟气象场的影响，分析区域平均值而不是单独的格点值，下文的讨论如无特别说明均针对沈阳城市20km×20km 区域(图 1b)平均值.0600LST 和1400LST 位温廓线如图 3所示，0600LST 出现接地逆温层，MYJ方案模拟低层气温低于YSU 方案，导致逆温强度更强，500 m 以下逆温强度:MYJ 1.5K/100m,YSU0.9K/100m.逆温层结抑制污染物的垂直扩散，所以MYJ模拟的气象场更容易导致污染物在低层累积.1400LST 地表受太阳辐射升温，加热低层大气，接地逆温消失，低层出现中性或弱不稳定层结，500 m 高度以上仍覆盖着悬浮逆温层.

图 3

Fig. 3

图 3 位温廓线(a)0600LST 和(b)1400LST Fig. 3 Potential temperature profile at (a) 0600 LST and (b) 1400 LST

为进一步分析温度层结日变化情况，图 4和图 5分别给出MYJ、YSU 模拟的位温和两种方案模拟差值(YSU-MYJ)的高度-时间剖面.MYJ方案中接地逆温维持时间较长，中性层结只在1400LST 左右维持两个小时.YSU 方案中接地逆温从1200 LST 开始消失，一直到1800LST 才重新建立，期间6个小时维持中性层结，使得低层积累的污染物有更多的时间向高层大气扩散.从图 5 可以看到，YSU 方案相比MYJ方案的增温主要出现在夜间低层大气，在日出前后0700LST 达到最强3 ℃左右. 日出后随着边界层发展，增温的垂直范围逐渐升高，在1400LST 达到最高700m 左右.

图 4

Fig. 4

图 4 不同边界层方案模拟位温的高度-时间剖面(a)MYJ；(b)YSU. Fig. 4 Height-time section of the potential temperature simulated with different boundary layer schemes

图 5

Fig. 5

图 5 YSU 方案和MYJ方案模拟位温差值的高度-时间剖面 Fig. 5 Height-time section of the difference between the potential temperatures simulated with YSU and MYJ

低层大气存在两种加热机制:地表接收太阳辐射增温，通过湍流热通量加热低层大气；位温较高的高层大气通过卷挟作用进入边界层，混合加热低层大气.YSU 方案模拟低层气温高于MYJ方案，为分析这种差异主要受上述哪种加热机制作用，研究地表感热通量的日变化.如图 6 所示，感热通量在0900LST 至1500LST 为正值，即热量从地表向大气输送，并在1200LST 达到最大100 W·m^-2 左右.期间MYJ方案模拟值都大于等于YSU 方案，同时2m 高度温度YSU 模拟一直高于MYJ，Hu等人^[15]的模拟得到类似结果，因此可以推断这种温度差异不是由地表加热造成的，而是受高层大气热量向下输送的影响.夜间地表辐射降温，温度低于低层大气，热量从大气向地表输送，即感热通量为负值，期间YSU 模拟的绝对值均大于MYJ，即更多的热量从低层大气输送到地表.夜间YSU 模拟低层气温一直高于MYJ，而且温差不断扩大，直到日出前后0700LST 达到最大(图 2和图 5)，相应的感热通量差值也在此时达到最大.YSU 方案相比MYJ方案，低层大气维持较大的向下感热通量，同时两方案间的模拟温差不断扩大，可见YSU 中低层大气通过高层大气的卷挟混合作用获得更多的热量.

图 6

Fig. 6

图 6 感热通量日变化 Fig. 6 Diurnal variation of the sensible heat flux

日出后0900LST 到1100LST 的位温变化廓线如图 7所示.随着太阳辐射增强，地表对低层大气的加热作用加强，MYJ中受较大的向上感热通量作用(图 6)，低层大气增温快于YSU.YSU 中垂直湍流混合更加有效，所以变温分布更均匀.注意到在1000m 高度附近存在负变温区，这是由于边界层顶的卷挟混合作用降低了高层位温，由于MYJ低层温度较低，导致高层负变温更大.Hu等人^[15]模拟夏季的情况发现在1000 m 高度以上存在负变温区，但YSU 变温比MYJ大.综上所述，YSU 方案模拟低层气温高于MYJ方案的原因是:YSU 边界层顶卷挟作用和边界层内垂直混合更强，热量从高层大气往下输送，加热低层大气，夜间由于没有地表加热作用，两模式间的模拟温差不断扩大，直到日出前达到最大.

图 7

Fig. 7

图 7 0900LST 至1100LST 位温变化廓线 Fig. 7 Profile of potential temperature change from 0900 LST to 1100 LST

水平风速决定了污染物的输送效率，静风小风条件往往造成污染物在排放源附近累积，同时风速较小时，湍流运动的量级更接近平均运动，加大准确预报风场的难度.为检验WRF 对不同风速的预报能力，图 8给出模拟和观测风速的比值随观测风速的分布.可以看到风速较小时，结果分布较散乱，随着风速增大，模拟观测比值逐渐收敛到y=1 线附近，图中两条虚线分别表示y=2和y=0.5，即模拟观测比值2倍线，落在2倍线内数据的比例(FAC2) 可以用来衡量模式的模拟能力.

图 8

Fig. 8

图 8 风速模拟和观测比值随观测风速分布 Fig. 8 Distribution of the ratio between simulated and observed wind speed as a function of observed wind speed

图 9给出FAC2随观测风速的变化，当风速小于4m·s^-1，FAC2 随着风速增大近似线性增长，对4m·s^-1以上风速，两种边界层方案的FAC2 趋于稳定，大部分都在90% 以上，对2 m·s^-1以下风速WRF模拟能力较差，FAC2 低于60%，1 m·s^-1以下的模拟甚至降到30%.空气污染多发生在稳定天气形势的小风条件下，WRF 对低风速的模拟能力还有待提高，不同边界层方案对风速模拟结果影响不大.

图 9

Fig. 9

图 9 风速FAC2随观测风速分布 Fig. 9 Distribution of FAC2 for wind speed as a function of observed wind speed

通过以上分析可以看到WRF模拟低层风速偏大，对低风速的模拟能力较差，而污染多发生于低风速条件下，因此进一步分析观测和模拟的风速累积概率分布情况(图 10).在模拟时段内，天气形势较稳定，观测风速中低于2 m·s^-1的占近60%，模式低估了小风出现的频率，只有30%的模拟风速低于2m·s^-1；4 m·s^-1以上风速在观测中占10%，在模拟中占30%.WRF模拟低估小风发生频率，高估大风发生频率，可能会使与之耦合的空气质量模式模拟的污染物扩散过快，导致浓度偏低.

图 10

Fig. 10

图 10 风速累积概率 Fig. 10 The cumulative probability of wind speed

风向是影响污染物扩散的重要因素，尤其在发生危险品泄漏事故时，决定了毒害区域的位置.表 2 给出观测和模拟的风向频率分布以及各风向对应的平均风速，模拟时段内观测结果表明北风出现频率最高，占21%，南风次之，西风和西北风频率最低； 北风和西北风平均风速最大2.9m·s^-1，东风风速最低1.5m·s^-1.两种边界层方案模拟的北风频率均为最高，但东风和东南风频率最低，较观测低七八个百分点.平均风速模拟结果偏大，但随风向的分布跟观测一致:北风风速最大，东风最小.

表 2(Table 2)

表 2 风向频率及对应平均风速 Table 2 The frequency of wind direction and the corresponding average wind speed

表 2 风向频率及对应平均风速 Table 2 The frequency of wind direction and the corresponding average wind speed

边界层高度(PBLH)是分析边界层结构的重要物理量，决定了可供污染物扩散稀释的潜在空气体积.PBLH 高时，污染物可以在更大的空间内扩散稀释，从而降低浓度.因此，PBLH 也是影响空气质量的重要指标，在空气质量数值模式中，需要由观测数据或气象模式提供PBLH，用于湍流扩散参数化^[41].

MYJ方案和YSU 方案用不同的方法计算PBLH:MYJ中PBLH 定义为湍流动能小于临界值0.1m²·s^-2 的最低高度(以下简称TKE 方法)；YSU 中定义为整体Richardson数小于临界值0 的最低高度(以下简称Ri-0方法).通过Richardson数计算PBLH 的方法运用较广泛，但在有些研究工作中临界值取0.25(以下简称Ri-0.25 方法)^[42-43].整体Richardson数定义如下:

(5)

其中θ_vk表示第k层和最低层之间的大气平均位温，ΔU_k，ΔV_k平均风速在第k层和最低层之间的差值，Δθ_vk位温差值，g重力加速度，Δz_k垂直距离.除上述TKE 和Ri方法外，Srivastava等^[44]根据位温廓线将PBLH 定义为垂直减温率大于等于2 K/km 的最低高度(以下简称PT 方法).有必要考察各种计算方法对PBLH 计算结果的影响.

表 3给出基于不同边界层方案模拟数据，采用不同方法计算得到的PBLH 情况，其中标准差表示PBLH 的日变化强度.可以看到，根据相同的PBLH 计算方法(Ri-0.25，PT)，YSU 方案模拟的PBLH 平均高度、日最大高度、标准差均大于MYJ方案. 垂直混合和边界层顶卷挟作用越强，越有利于边界层垂直发展，因此近一步证实YSU 方案垂直混合和卷挟作用比MYJ方案强.比较不同的计算方法可以看到，在MYJ方案中TKE 方法计算的PBLH 最大，PT 方法最小，两者平均高度相差接近100m,日最大高度均值相差超过300m.在YSU 方案中未计算湍流动能，所以没有采用TKE 方法.YSU 的三种计算方法中，Ri-0.25 和PT 的平均高度较接近，但Ri-0.25的最大高度大于PT，两者相差超过150m,Ri-0.25 的标准差也大于PT，可见中午采用Ri- 0.25方法计算的PBLH 高于PT，导致PBLH 日变化增大.这可能是由于Ri-0.25 综合考虑了热力驱动和动力剪切作用对边界层发展的贡献，而PT 中只考虑热力因素.比较YSU 方案中Ri-0和Ri-0.25 方法的计算结果发现，临界Richardson数从0提高到0.25，对应的PBLH 也变大，平均高度增加50m 左右，最大高度增加近200m,同时日变化也更加强烈.可见通过Richardson数计算PBLH 的方法对临界值的选取较敏感，需要分析更多的观测数据确定合适的临界值.比较边界层方案自身输出的PBLH，MYJ方案(基于TKE)大于YSU 方案(基于Ri-0)，和Borge等人^[19]的模拟结果一致.

表 3(Table 3)

表 3 不同边界层方案不同计算方法得到的边界层高度 Table 3 The boundary layer height derived from different boundary layer schemes with different methods

表 3 不同边界层方案不同计算方法得到的边界层高度 Table 3 The boundary layer height derived from different boundary layer schemes with different methods

本文运用中尺度气象模式WRF 模拟沈阳地区冬季边界层结构，分析两种边界层方案MYJ 和YSU 对模拟结果的影响，深入研究跟污染物扩散密切相关的温度层结、低层风场、边界层高度等物理量，得到以下初步结论:

(1) YSU 方案中边界层顶卷挟和边界层内混合作用强于MYJ方案，热量向下输送更加有效，导致模拟低层气温高于MYJ，最大温差出现在日出前后；夜间YSU 模拟的逆温强度低于MYJ，逆温维持时间比MYJ短4小时，有利于污染物垂直扩散.

(2) 风速风向对边界层方案的敏感性不如温度明显.WRF对2m·s^-1以下风速的模拟较差，模拟风速偏高，低估低风速出现的频率，高估高风速出现的频率，这种偏差会导致空气质量模式模拟浓度偏低.研究时段内，沈阳区域以北风为主，对应风速较大，东风对应风速较小，WRF 模拟的风向风速对应关系和观测一致.

(3) YSU 方案模拟的边界层高度大于MYJ方案，更有利于污染物的扩散.不同边界层高度计算方法中，通过湍流动能计算的高度最高，整体Richardson数次之，通过位温计算的高度最低.其中Richardson数方法中的临界值选取对结果影响较大，需进一步通过观测验证.

本项研究揭示的初步结果，可以在气象模式和空气质量模式的耦合研究中得到更加深入的解释应用，为城市环境健康和谐发展提供科学决策依据.