北京冬季城市人为水汽排放的定量估算研究

徐一平 苗世光 赵德龙 丁德平

1 中国气象科学研究院，北京 1000812 北京城市气象研究院，北京 1000893 北京市人工影响天气办公室，北京 1000894 云降水物理研究和云水资源开发北京市重点实验室，北京 1000895 中国气象局华北云降水野外科学试验基地，北京 101200

1 引言

城市化是人类社会现代化发展的必然趋势，它使得土地利用和土地覆盖都发生了很大的变化，人为热排放、城市气溶胶和人为温室气体排放也大大增加，对天气、气候和环境均产生了不同程度的影响，天气、气候和环境反过来也会对城市的运行和居民的生产生活产生影响，城市化程度体现了一个国家的经济发达程度，是工业化水平高低的重要指标（苗世光等, 2014b）。截至2019 年底，我国城市化率已达到60.6%，城镇人口突破8.4×109（中华人民共和国国家统计局, 2020）。

目前国内大部分的研究主要集中在城市热岛、城市大气环境上（汪光焘等, 2005; 寿亦萱和张大林, 2012; Deosthali, 2000），而对于城市人为水汽研究的相对较少的。国外在这方面已经开展了一些相关研究。Moriwaki et al.（2008）使用地理信息系统和能源消耗数据库的数据，估算了东京人为水汽排放量和人为热的时空变化，发现夏季东京市中心的人为水汽最大值超过500 W/m-2。冬季东京居民区的人为热排放量大于夏季商业区的人为热排放量，并且人为热和水汽，使气温升高了约0.5°C，水汽混合比增加0.1 g/kg（Senoo et al., 2004）。在美国进行的飞机观测实验中，在城市下风方向观察到较高的水汽摩尔分数，水汽的排放量的估计值在1.6×104和1.7×105kg/s之间，占城市研究点大气边界层总的水汽平流的8.4%（Salmon et al., 2017）。

人为水汽的排放也会对城市的天气、气候和环境等带来多方面的影响。天然气燃烧产生的水汽会使空气饱和，在低温的共同作用下，产生了埃德蒙顿的低温雾（Rage, 1972）。城市的过量水汽与城市热岛之间也存在着正相关，但是城市过量水汽对城市热岛的增强作用较弱，并且存在相当大的负反馈（Holmer and Eliasson, 1999）。通常城市热岛最强时，经常观察到城市过量水汽（Kuttler et al.,2007）。另一方面，来自城市的过量水汽也会导致城区下风向的降水增加（Mölders and Olson, 2004）。由此可以看出，城市水汽的排放对于城市的天气、环境等也有不可低估的影响。

作为我国的首都和政治文化中心，北京市人口规模大，经济活动发达，属于特大城市，近几十年来，北京的城市建设十分迅速，经历了高速发展期，城市结构和建筑高度等都发生巨变（徐阳阳等,2009）。随着城市规模迅速扩大和城市化步伐的加快，北京的气象条件发生了一系列的变化。因城市效应，北京冬季表现为弱湿岛，城区与郊区比湿日变化有明显差异（窦晶晶等, 2014）。北京城区冬季寒冷干燥，波恩比较大，白天储热在净辐射中所占比例最大，人为热中潜热和感热表现为同等重要，多云时，冬季夜间的感热略大于潜热（苗世光等,2012）。在城市区域，人为潜热排放与人类活动的规律基本吻合，比如空调系统、居民生活排放和汽车尾气排放等均会导致人为潜热的大量排放（苗世光等, 2014a）。进一步研究发现，城市对于北京地区的湿度影响显著（Liu et al, 2009）。总体上来说，城市在白天会比周围地区更加干燥（Arnfield, 2003）。

根据北京市统计局资料显示，2019年北京市常住人口2153.6×104，常住人口密度达到1312 km-2（北京市统计局, 2020）。人口密度大，人类活动强烈，人为水汽排放量巨大，十分具有研究意义。本文利用飞机观测数据量化了北京整个城市内水汽摩尔分数的增强情况，并讨论了可能的来源，将为数值模式中人为潜热的设置提供参考。

2 研究方法

本研究采用筛选出的2016年12月16日、17日、18日、19日符合质量平衡飞行的飞机观测数据，计算了城市过量水汽通量、城市过量水汽排放率和城市水汽总排放率。

2.1 航线位置

北京市气象局人工影响天气办公室的飞行观测实验，每次均从北京北部郊区出发，环绕北京中心城区飞行观测，并飞行多圈，最后返回北部郊区，完成一次飞行观测实验，每次飞行的航线基本类似。

2.2 机载平台和机载仪器

2.2.1 机载平台和机载仪器的具体型号

本文观测仪器的机载平台为北京市人工影响天气办公室空中国王3587飞机。本文使用的数据来自机载气象探测系统AIMMS-20（Aircraft-Integrated Meteorological Measurement System-20）。该系统可用于飞行时的温度、湿度、风的三维空间状况。该仪器包括空气数据探头，通过压差和逆流温度/湿度测量室测量气流和空气温度湿度；惯性测量单元测量整个航迹中飞机高度数据；GPS测量飞机高度和机身倾斜角度。

2.2.2 机载仪器的指标

风速精度水平分量为0.50 m/s，垂直分量为0.75 m/s；温度精度为0.3°C，分辨率为0.01°C；相对湿度精度为2%，分辨率为0.1%；工作温度为-20～50°C。

2.3 数据

本文使用数据是北京市气象局人工影响天气办公室3830飞机观测资料，筛选出温度、相对湿度、气压、位置（经纬度）、飞行高度、东西向和南北向飞行速度等数据。因组织一次专门飞机观测实验的成本较高，本文研究所选用的数据是从北京市气象局人工影响天气办公室的日常作业观测的3830飞机观测数据中筛选出来的，因此存在数据针对性不强，很多数据不满足计算条件而无法使用的情况，但是通过一定的筛选条件，依然能够从广泛的数据中可以找到可使用的部分，为本次研究提供数据基础。

2.4 质量平衡飞行

筛选出的飞机航线所用时间均为数分钟，以尽量减少飞行期间大气边界层的增长。飞机在北京市的上风和下风向飞行，城市过量水汽排放率根据中心城区相对于背景的摩尔分数的增加、相对于飞行轨迹的风速垂直分量和当时的边界层高度计算。图1、图2 以2016年12月16日的一次飞行为例，显示了飞机观测航线和飞行高度变化的大致情况，所有日期飞行观测的航线大致类似，并且在该次飞行中观测到了过量水汽排放信号，即人为水汽排放。

图1 2016年12月16日的飞机观测航线其中一圈。灰色区域表示城市区域，红色曲线表示航线，黑色环线为北京二至五环路环线Fig. 1 One circle of the aircraft observation route on 16 December 2016; the gray area represents the urban area, the red curve represents the route,and the Beijing Second to Fifth Ring Road are drawn

图2 2016年12月16日的一次观测飞行中飞行高度随飞行时间变化Fig. 2 Change of flight altitude with flight time in a flight on 16 December 2016

2.5 背景浓度的确定

该研究的目的是确定城市过量水汽排放率。背景浓度可作为确定中心城区水汽摩尔分数增量的参考。本次研究中采用中心城区周边区域作为背景，周边区域人口聚集度较小，相对中心城区，人为潜热排放量也较小。具体计算的过程中，即根据位于城区下风向航线两侧边缘的数据确定背景水汽摩尔分数，背景浓度与水平方向飞行距离线性相关。

2.6 天气形势

根据筛选出的航线的飞行时间，从北京市气象台获取了2016年12月16～19日每天08:00（北京时间，下同）海平面、850 hPa、700 hPa、500 hPa的天气图。从海平面气压场等压线分布图可见，12月16日北京地区处于弱高压前部的均压场中，弱高压的中心位于北京的西部；12月17日高压系统东移入海，北京处于两个高压之间的低压辐合区控制；12月18日北京持续处于两个高压之间；12月19日北京上游高压逐渐东移南下，北京处于高压前部的弱气压场中。大气层结垂直方向上，16～19日，500、700、850 hPa高度上，北京均位于槽后，为下沉气流。总体看来，2016年12月16～19日，北京受高压控制，无降水，水平方向风速较小，无明显暖湿气流，天气稳定为主。

2.7 城市过量水汽通量、城市过量水汽排放率、城市水汽总排放率的定义和计算方法

城市过量水汽通量是指单位时间内，水平方向上流经单位面积城市比周围乡村多的水汽量，单位用mol m-2s-1表示；城市过量水汽排放率指的是单位时间内，城市比周围乡村多排放的水汽质量，单位用kg/s表示；城市水汽总排放率指的是单位时间内，整个城市排放的水汽质量，单位用kg/s表示。对于城市过量水汽通量FUrban,i，计算公式如下：

其中，i表示每一个下风向数据点的标号，U⊥,i表示垂直于航线的水平风分量，cdw,i为下风向水汽摩尔分数，cbg,i为背景水汽摩尔分数，使用理想气体定律，可将水汽摩尔分数转化成水汽摩尔密度。相同气压和温度下，下风向水汽摩尔密度减去相应的背景水汽摩尔密度，是城市过量水汽摩尔密度，城市过量水汽摩尔密度乘以风速相对于航线的垂直分量就可以得到城市过量水汽通量FUrban,i，同理，城市水汽通量FTotal,i计算公式如下：

为了计算城市范围内的水汽排放率，需要知道城市过量水汽通量矩阵，并做积分

其中，H为大气边界层高度，x为水平飞行距离的一半，M是指城市过量水汽通量矩阵或者城市水汽总通量矩阵，进一步比较可以确定城市过量水汽排放率（EUrban） 与城市水汽总排放率（ETotal）的比值。

3 计算过程

3.1 饱和水汽压计算

飞机观测的气压范围大致在640～1020 hPa，仍以0°C作为水的凝固点，饱和水汽压es的计算，温度在0°C以上使用水面饱和水汽压公式：

温度在0°C以下使用冰面饱和水汽压公式：其中，t为摄氏温度，T为热力学温度。

3.2 水汽摩尔分数的计算

其中，χV为水汽摩尔分数，e为水汽压，p为气压，UW为 相对湿度；es(T)为饱和水汽压，由公式（4）、公式（5）可得；χVS为饱和水汽摩尔分数。根据公式（6）和公式（7）则

3.3 将水汽摩尔分数转化为摩尔密度

利用水汽的理想气体状态方程：

其中，n表示水汽的物质的量，V表示气体体积，可得水汽摩尔密度，R=8.314 J mol-1K-1。

3.4 数据的筛选

北京供暖时间为11月15日左右，飞机飞行日期为2016年12月16日、17日、18日、19日，北京地区已完全进入冬季，之所以选择冬季的数据，是因为冬季的蒸发速率和饱和水汽压低于夏季，并且冬季人为水汽排放量大、背景水汽排放量小，易于估算城市人为水汽排放通量和排放率。

根据公式（1）～（3）判断，飞机观测数据的筛选条件主要包括以下6点：（1）飞机观测过程中，风向需指向城市外侧，即使得城市位于下风向； （2）飞机观测飞行高度需要在探测时大气边界层高度以下；（3）飞行路径需要围绕北京市区，郊区的飞行路径属于无效路径；（4）飞行路径的侧边界尽量远离城区；（5）飞机飞行路径可根据经纬度信息进行判断，当飞行轨迹不为大致南北、东西向时，需要考虑飞行轨迹与正北向的夹角；（6）数据需无缺测。利用以上的这些筛选条件，本研究对2016年12月16日、12月17日、12月18日、12月19日的飞机探测数据进行了筛选，最后得到了可用数据。

3.5 大气边界层高度的计算

在本研究中边界层高度的确定采用位温梯度法，位温梯度法计算边界层高度的依据是：近地层以上 （不稳定边界层时为混合层以上），位温梯度最大的高度定义为边界层高度。由于没有连续的位温资料，此处采用中央差分：

其中，n=2,3,4,···,θn+1表示n+1点处的位温，θn-1表示n-1处的位温， Δz表示n+1点与n-1点之间的距离，根据飞机探空曲线可算出边界层高度。之所以采用1次观测代表大气边界层高度，是因为在没有大尺度天气系统影响下，混合层盖顶逆温位置变化相对比较缓慢，1天中1次观测的盖顶逆温底高度可以代表这1天的情况，混合层顶覆盖逆温的位置就是大气边界层高度（袁仁民等, 2003）。

4 结果分析

2016年12月16日位温廓线和可用航线选择的结果如图3所示，利用上升阶段飞机观测到的气压、温度、高度数据，可绘制位温廓线，通过位温梯度法，920 m处位温梯度最大，为0.022 K/m，确定的大气边界层高度为920 m，并通过6点筛选条件确定了07:30:57至07:33:30自东向西、07:40:15至07:42:45自南向北的可用航线。

图3 （a）2016年12月16日北京的位温廓线；（b）本文选取的2016年12月16日07:30:57至07:33:30（红色）、07:40:15至07:42:45 （绿色）可用飞行路径Fig. 3 (a) Potential temperature profile on 16 December 2016 in Beijing; (b) the red and green are the selected available flight path maps during 0730:57–0733:30 BT (Beijing time) (red) and 0740:15–0742:45 BT (green) on 16 December 2016, respectively

12月17日位温廓线和可用航线选择的结果如图4所示，同样利用上升阶段飞机观测到的气压、温度、高度数据，通过计算可绘制位温廓线，通过位温梯度法，888 m处位温梯度最大，为0.053 K/m，确定的大气边界层高度为888 m，并确定了06:39:33至06:43:00有一定倾角的自北向南、07:03:38至07:06:16和07:06:58至07:10:26两条自西向东的可用航线。

图4 （a）2016年12月17日北京的位温廓线；（b）本文选取的2016年12月17日06:39:33至06:43:00（红色）、07:03:38至07:06:16（ 绿色）、07:06:58至07:10:26（黄色）可用飞行路径Fig. 4 (a) Potential temperature profile on 17 December 2016 in Beijing; (b) the red, green, and yellow are selected during 0639:33–0643:00 BT,0703:38–0706:16 BT, and 0706:58–0710:26 BT on 17 December 2016, respectively

12月18日位温廓线和可用航线选择的结果如图5所示，类似的，利用上升阶段飞机观测到的气压、温度、高度数据，绘制位温廓线，570 m处位温梯度最大，为0.031 K/m，确定的大气边界层高度为570 m，并确定了05:33:56至05:36:05较小倾角自北向南、05:37:04至05:39:57自西向东、05:40:32至05:43:24自南向北的可用航线。

图5 （a）2016年12月18日北京的位温廓线；（b）本文选取的2016年12月18日05:33:56至05:36:05（红色）、05:37:04至05:39:57 （绿色）、05:40:32至05:43:24（橙色）可用飞行路径Fig. 5 (a) Potential temperature profile on 18 December 2016 in Beijing; (b) the red, green, and orange are selected available flight path diagrams during 0533:56–0536:05 BT, 0537:04–0539:57 BT, and 0540:32–0543:24 BT on 18 December 2016, respectively

12月19日位温廓线和可用航线选择的结果如图6所示，通过位温廓线，617 m处位温梯度最大，为0.023 K/m，确定的大气边界层高度为617 m，并确定了06:02:49至06:05:38和06:16:42至06:19:28自西向东，06:03:31至06:09:17自南向北的可用航线。

图6 （a）2016年12月19日北京的位温廓线；（b）本文选取的2016年12月19日06:02:49至06:05:38（黄色）、06:03:31至06:09:17 （绿色）、06:16:42至06:19:28（红色）可用飞行路径Fig. 6 (a) Potential temperature profile on 19 December 2016 in Beijing; (b) the yellow, green, and red are the selected available flight path diagrams during 0602:49–0605:38 BT, 0603:31–0609:17 BT, and 0616:42–0619:28 BT on 19 December 2016, respectively

在确定了可以使用的观测数据之后，利用公式 （1）～（3），计算可得城市过量水汽通量（FUrban）、城市过量水汽排放率（EUrban）和城市水汽总排放率（ETotal）、城市过量水汽排放率（EUrban）占城市水汽总排放率（ETotal）的比例，如表1所示。

表1 飞行时间、城市过量水汽通量（FUrban）、城市过量水汽排放率（EUrban）和城市水汽总排放率（ETotal）及其比例Table 1 Flight time, urban anthropogenic water vapor flux (FUrban), urban anthropogenic water vapor emission rate (EUrban),and total urban water vapor emission rate (ETotal) and their proportions

在北京附近连续4天的飞机探测中选取的11次可用的飞行路径中，有5次（2016年12月16日07:30:57至07:33:30、2016年12月18日05:33:56至05:36:05、2016年12月18日05:40:32至05:43:24、2016年12月19日06:06:31至06:09:17、2016年12月19日06:16:42至06:19:28）观察到中心城区的城市水汽信号升高，可认为这些时间内中心城区相对于周边区域存在水汽的过量排放。在筛选出的中心城区水汽过量排放的实例中，表现为中心城区的水汽摩尔分数超过周边地区，过量水汽排放的占比在0.64%～13.12%。过量水汽排放的信号大小每时每刻都是在变化的，在有的飞行探测试验中，边界层的垂直混合程度也会产生影响。

从表格中可以看到在过量水汽排放的计算中，也会出现出一些负值，说明随着北京周边城市的快速发展，被作为背景浓度考虑的四环到五环之间的水汽排放也逐渐增多，甚至在某些时候已经超过中心城区的水汽排放量。可能原因是随着经济的发展，北京四环与五环之间成为新的经济中心，城市建设的快速发展与人口的大量集聚，导致人为潜热的排放大量增加。对于2016年12月16日07:40:15至07:42:45飞行中计算出的城市过量水汽通量和城市过量水汽排放率的负值，主要原因是其作为背景浓度的一个取值点位于北京的望京区域，该区域商业和写作楼密集，存在较大的人为水汽排放；另一个点位于北京经济技术开发区附近，产业密集，人口众多，同样存在较大的人为水汽排放，因此使得计算过程中确定的背景浓度相对较高，导致了城市过量水汽通量和排放率负值的出现。对于2016年12月17日06:39:33至06:43:00飞行中计算出的城市过量水汽通量和城市过量水汽排放率的负值，主要原因是作为背景浓度的一个取值点位于中关村科技园附近，经济和科技发达，人员和商业密集，人为水汽排放量巨大；另一个背景浓度取值点位于永定河附近，有众多旅游景点、企业和居民区，同样有较大的人为水汽排放量，同16日的飞行类似，计算中确定的背景浓度相对较大，导致了负值的出现。其他负值的分析类似。

本文根据城区下风向航线两个侧边缘的数据确定背景水汽摩尔分数，若航线侧边缘处城市化水平较高，人为水汽排放量大于中心城区，就会导致负值的出现。

如果有专门的飞行试验，使得飞机航线的两个侧边缘远离城区和大型水体，更好地代表乡村，则不会出现负值。

将所有观测到水汽过量信号的航线及所对应的城市过量水汽排放率占城市水汽总排放率的比例绘制在一张图上，如图7所示，在观测到水汽过量信号的航线中，北京城市北部地区的城市过量水汽排放率占城市水汽总排放率的比例显著大于北京城市南部地区，而东西两侧差异较小，在城市水汽总排放率相差不大的前提下，表明北京城市北部地区的人为潜热排放显著大于北京城市南部地区。

图7 2016年12月16日07:30:57至07:33:30（红色）、18日05:33:56至05:36:05（黄色）、18日05:40:32至05:43:24（绿色）、19日06:06:31至06:09:17（蓝色）、19日06:16:42至06:19:28（橙色）测得城市过量水汽信号的航线及由此计算出的城市过量水汽排放率占城市水汽总排放率的比例（EUrban/ETotal）Fig. 7 Route of the measured urban water vapor excess signal during 0730:57-0733:30 BT (red) on 16 December 2016, 0533:56-0536:05 BT(yellow), 0540:32-0543:24 BT (green) on 18 December 2016, 0606:31-0609:17 BT (blue), 0616:42-0619:28 BT (orange) on 19 December 2016 and the ratio of the calculated urban water vapor excess emission rate to the total urban water vapor emission rate (EUrban/ETotal)

对于北京城市北部地区人为潜热大于南部的现象，主要从天气背景和城市内部结构两个方面分析。天气背景方面，如果飞行当日北京相对湿度较大，人为水汽占比小，计算会出现一定的误差。城市内部结构方面，这可能跟北京城市北部城市化程度高、人口密度较大有关。

5 结论与讨论

本文利用飞机观测资料对北京冬季人为水汽排放进行了定量估算，主要计算了城市过量水汽通量、城市过量水汽排放率、城市水汽总排放率以及城市过量水汽排放率占城市水汽总排放率的比例，得到如下的主要结论：

（1）在北京城区下风方向进行的飞机观测中多次观测到较高的水汽摩尔分数，当观测到城市水汽过量信号时，水汽排放率的估计值在57.64到9538.87 kg/s之间，占城市大气边界层水汽总平流 （背景+城市过量）的0.64%到13.12%。

（2）城市下风向观察到的水汽摩尔分数过量的原因，可能有以下3个方面，第一是人为的水汽来源，如燃烧和蒸发冷却，近年来北京市推行的煤改气政策可能起到了增加中心城区人为水汽排放的作用；第二是局地水体的贡献，如玉渊潭等城区水体；第三是中心城区和周边地区融雪率和蒸发率的差异，在本研究中，这个差异较小。

本研究观测到的城市过量水汽排放发生在一年中蒸发速率比较低的冬季，而且，由于较低的气温，冬季的饱和水汽压低于夏季，相对于噪声和自然变率，城市水汽过量信号更加明显。使用机载平台在北京城市下风方向对过量的水汽摩尔分数的研究，显示的是全市范围内的水汽浓度，而不是局部的或者相对来说较小尺度的观测结果。

由于飞行航线不是专门为本研究设计，因此存在可用数据量少、筛选困难等情况。如果要进一步对该问题进行研究，可设计专门的飞行试验，进一步扩大航线范围，使得背景浓度更能代表受城市影响较小的区域；在计算水汽排放率时也可以增加垂直方向上高度的变化等。另一方面，本文使用的飞机观测数据是从北京市气象局人工影响天气办公室日常飞机探测的数据中筛选出来的，没有设计和进行专门的飞行观测试验，成本较低，可操作性强，为在国内不同城市开展城市尺度人为潜热的计算提供了新的思路和方法，最终的计算结果也为城市气象环境数值模拟中人为潜热的设置提供了一个参考。