基于WRF-CMAQ的中国典型湖泊大气氮沉降特征

樊洋成,刘 萍

(1.上海交通大学 中英国际低碳学院,上海 201306; 2.上海交通大学 环境科学与工程学院,上海 200240)

0 引 言

近2个世纪以来,由于全球规模的工业化和氮肥、化石燃料的大量使用,人为源活性氮(N)排放量大幅上升[1]。过去30 a里,中国每年的活性氮沉降总量增加了约25%,从11.1 kgN/hm2增加到了15.3 kgN/hm2[2],导致内陆水体中的氮素浓度每年约增加76.6～93.9 μg/L。大量的大气氮沉降会加剧一系列环境问题,包括湖泊氮污染。一旦湖泊水环境中的氮浓度超过阈值,就会发生水体酸化以及富营养化现象,进而引发生物多样性的下降以及水生态环境恶化。国外有研究表明,加拿大苏必利尔湖的氮素浓度近50 a持续上升,大气氮沉降是主要原因之一[3];国内有研究指出,2019年安徽巢湖湖面大气总氮年沉降量达到661.16×103kg,占巢湖主要河流入湖负荷的48.41%。大气氮沉降是当地农田生态系统和巢湖水体氮输入的重要来源[4]。

随我国经济的快速发展,越来越多的活性氮由工业源与农业源排放到大气中,全国年均总氮沉降通量从20世纪60年代到近10 a增长了约59%[5]。华东地区的无机氮年沉降通量约为15 ～ 50 kgN/hm2,远高于我国西部的平均水平(约2 kgN/hm2)[6]。国内针对不同湖泊的已有研究表明,滇池[7]、洞庭湖[8-9]、太湖[10-14]、青海湖[15]、东湖[16]、洱海[17-18]、乌梁素海[19]流域的大气氮沉降与当地的氮肥施用、畜禽养殖、化石燃料消耗以及生物质燃烧等人类活动有不同程度的关联性,然而大气氮沉降在我国不同湖泊的时空分布特征及对湖泊水体氮污染贡献程度的研究还相对较少。目前的研究多针对某一区域,实地设点采样,通过对降尘与雨水样品进行化学定量分析的方式对氮沉降通量进行测量[4,20]。对于较大区域内的沉降情况,由于难以获取高精度和连续的观测数据,多采用大气传输模型进行数值模拟,常用于研究氮沉降的模型包括CMAQ、EMEP Unified、FRAME等[21]。本研究采用中尺度天气预报系统及区域多尺度空气质量模式(WRF-CMAQ)[22-23]对2018年中国大气氮沉降通量进行了模拟。WRF(The Weather Research and Forecasting)模式由美国国家大气研究中心等多方开发,广泛应用于天气研究与业务预报[24],本研究使用其为CMAQ模式提供气象场数据。CMAQ(Community Multiscale Air Quality system)模式由美国环境保护署(Environmental Protection Agency,EPA)研发,多为大气污染控制策略提供决策支撑,国内外已有较多研究使用CMAQ模式对大气氮、硫等污染物元素沉降通量进行模拟研究,模拟结果与实测值普遍存在一定差异,但较长时间尺度内的变化趋势基本吻合[25-27]。

本研究聚焦我国不同地区的4个典型湖泊(青海湖、乌梁素海、东湖、太湖),估算了4个湖泊区域不同形态氮的湿沉降、干沉降以及总沉降氮通量,分析了各湖泊沉降通量的时空变化特征和大气氮沉降对于各湖泊总氮负荷的贡献比例,探讨了影响大气氮沉降通量的主要因子,以期为通过制定污染控制措施减少大气氮沉降及湖泊水体氮污染提供数据支持和科学依据。

1 研究方法

1.1 研究区域

模拟区域以(31°N,120°E)为中心坐标,覆盖我国绝大部分区域,4个湖泊位置详情如下:青海湖(37°N,100°E)面积4 625.6 km2,位于地广人稀的青藏高原,是我国面积最大的内陆咸水湖[28];乌梁素海(41°N,109°E)位于华北内蒙古干旱草原与荒漠地区,面积<300 km2[29];东湖(31°N,114°E)是我国最大的城中湖,位于长江中下游平原的大型城市中心,面积约32 km2[30];太湖(31°N,120°E)位于长三角发达地区,邻近多个工农业相对发达的城市,面积约2 445 km2[31]。

1.2 模式参数配置

本研究使用WRF v4.2.1模拟气象场数据,中心坐标点为(31°N,120°E),垂直方向16层,水平分辨率12 km×12 km,模式顶大气压为50 hPa,网格数为512×512,初始场和边界场使用美国国家环境预报中心(National Centers for Environmental Prediction, NCEP)的再分析资料FNL(Final Operational Global Analysis)[32],模拟时间为2018年1—12月,模拟区域为4个典型湖泊,物理参数化方案见表1。2018年4个典型湖泊所在城市(青海湖、乌梁素海、东湖和太湖分别对应西宁市、巴彦淖尔市、武汉市和常州市)均未发生极端天气事件,2018年降水总量分别为518.9、196.4、1 110.6、1 204.2 mm,2001—2020年平均降水总量分别为441.4、240、1 269、1 091.4 mm[38],四地2018年降水总量与2001—2020年平均降水总量相差均<20%,气象因素具有代表性。

表1 WRF模式物理参数化方案Table 1 Physical parameterization schemes used in the WRF model

本研究使用CMAQ v5.3.1模拟大气污染物的排放、生成、转化、传输、沉降和去除等物理和化学过程。气相化学反应机制选用SAPRC-07[39],气溶胶机制选用AERO7[40](Aerosol Module Version 7)。WRF模式模拟的气象场,通过气象化学界面处理器(Meteorology Chemistry Interface Processor,MCIP)处理为CMAQ所需的网格化数据。模拟所用排放清单来自清华大学建立的中国多尺度排放清单模型(Multi-resolution Emission Inventory for China,MEIC)[41-42]与亚洲人为源排放清单MIX[43],其中MEIC清单为2017年数据,覆盖中国大陆区域,MIX清单为2011年数据,覆盖模拟区域内中国大陆以外的地区。CMAQ模式模拟区域与水平网格精度设置与WRF模式一致,垂直分为14层,层顶高度100 mb。初始场和边界场由Geos-Chem模型模拟的全球大气污染物浓度提供。

1.3 模拟效果评估方法

为评估模式模拟效果,将WRF模式模拟得到的站点的地面2 m温度(T2)、地面2 m相对湿度(RH2)和地面10 m风速(WSPD10)的小时均值数据与中国气象数据网[38]的相应观测数据进行比较,将由 CMAQ模式模拟得到的站点日最大臭氧(O3)8 h平均浓度、二氧化氮(NO2)1 h平均浓度、细颗粒物(PM2.5)24 h平均浓度与中国环境监测总站[44]的污染物浓度观测数据进行对比,通过统计分析评估模式模拟结果的准确性和可靠性。

本研究使用标准化平均偏差 (Normalized Mean Bias,NMB)、标准化平均误差(Normalized Mean Error,NME)以及皮尔森相关系数(R)来评估模式模拟效果。NMB、NME与R的计算公式分别为:

(1)

(2)

2 模拟结果验证

2.1 气象参数模拟效果评估

本研究对比了2 421个站点T2、RH2、WSPD10的模拟值和观测值,T2的NMB与NME分别为-4%与16%;RH2的NMB与NME分别为-5%与17%,WSPD10的NMB与NME分别为44%与70%,表明气象场模拟效果较好。将4个典型湖泊所在城市2018年的T2、RH2、WSPD10模拟小时均值与观测数据进行对比,评估结果如表2所示。

表2 典型湖泊区域WRF气象参数模拟效果评估Table 2 Evaluation of simulated meteorological parameters of WRF over typical lake areas

图1为2018年1、4、7、10月(分别代表春、夏、秋、冬4个季节)T2、RH2、WSPD10的日均模拟值与观测值的时间序列。4个城市T2观测和模拟值吻合得较好,相关系数R都在0.97以上,标准化平均偏差为-1%～25%;RH2的相关系数R都在0.80以上,标准化平均偏差为-7%～-3%。WSPD的R为0.40～0.70,NMB为-7%～23%,4个城市的WSPD10模拟值都较好地反映了观测值的变化趋势。武汉市与常州市均位于我国长江中下游平原地区,WSPD10的NMB分别仅为-7%和2%。西宁市与巴彦淖尔市的WSPD10模拟值略高于观测值,模拟偏差可能与我国西北与华北高原地区复杂的下垫面与地理条件有关。

图1 2018年1、4、7、10月4个城市T2、RH2、WSPD10模拟值与观测值日均值时序结果Fig.1 Time series of observed and simulated daily average 2-m temperature (T2), 2-m relative humidity (RH2), and 10-m wind speed (WSPD10) in four cities in January, April, July, and October in 2018

2.2 O3、PM2.5、NO2浓度模拟效果评估

本研究对比了全国31个省会国控站点最大O38 h平均浓度、NO21 h平均浓度、PM2.524 h平均浓度的模拟值与观测值,O3的NMB和NME分别为13%和34%,NO2的NMB和NME分别为-4%和58%,PM2.5的NMB和NME分别为-41%和52%,均符合生态环境部发布的《环境空气质量模型遴选工作指南》[45]中的对应评价标准(O3:-15%

图2 2018年1、4、7、10月O3、NO2与PM2.5浓度模拟值与观测值日均值时序结果Fig.2 Time series of simulated and observed daily O3, NO2, and PM2.5 concentrations in four cities in January, April, July, and October in 2018

表3 O3、PM2.5和NO2浓度模拟效果评估Table 3 Evaluation of simulated O3, NO2, and PM2.5 concentrations over typical lake areas

4个城市的最大O38 h平均浓度模拟值与观测值基本吻合,NMB为-6%～25%,相关系数R为0.50～0.80。NO21 h平均浓度NMB为-28%～44%,相关系数R为0.44～0.53。武汉市与常州市的PM2.5524 h平均浓度模拟值较好地复现了观测值,NMB分别为-1%与-32%,相关系数R分别为0.80与0.75。西宁市与巴彦淖尔市的PM2.524 h平均浓度模拟值低于观测值,NMB分别为-81%与-70%,可能是由于两地10 m风速的模拟值高于观测值,大气扩散作用被高估导致。此外,排放清单的不确定性会导致模拟结果产生一定程度的偏差[46-47],例如未考虑当地的沙尘天气以及交通工具、建设工地引发的扬尘等因素,也可能引起模式对西部城市PM2.524 h平均浓度的低估。本研究中WRF-CMAQ对于化学场的模拟效果和偏差水平与其他研究结果相近[48],表明模拟结果可以用于进一步分析大气氮沉降通量特征和评估大气氮沉降对湖泊氮沉降的影响。

3 结果与讨论

3.1 4个湖泊大气氮沉降通量特征分析

模拟区域内总氮(TN)、氧化态氮(OXN,主要包括NO2、NO、无机硝酸盐、无机亚硝酸盐、有机氧化态氮等)、还原态氮(REDN,主要包括NH3、无机铵盐、有机胺类等)的季度沉降通量空间分布见图3。结果表明,在华中地区与东南部沿海地区,氧化态氮与还原态氮在夏季的沉降通量都高于其他3个季节。夏季华北平原、华中地区还原态氮沉降通量尤为集中,季度累积量为8 ～16 kgN/hm2,与其他研究结果一致[8, 9, 49-50]。我国还原态氮排放绝大部分来自农业源畜禽养殖及氮肥施用,华北平原、华中、华东地区均为农业大区,夏季较为频繁的农业生产活动及前体物的挥发造成了上述区域较高的氮沉降通量,也会为位于上述区域的湖泊带来较高的氮负荷。氧化态氮沉降通量的空间分布显示出与城市关联较为紧密,春、夏、秋、冬四季中,长三角、珠三角城市群是氧化态氮沉降的主要区域,4个季节的累计氧化态氮沉降通量均超过5 kgN/hm2,由于大气氧化性增强,大气氧化态氮沉降对于临近城市湖泊的氮负荷贡献也会增加。

注:春季为3—5月份,夏季为6—8月份,秋季为9—11月份,冬季为12—翌年2月份。图3 2018年四季氮沉降通量空间分布Fig.3 Spatial distributions of simulated seasonal nitrogen depositions in four seasons in 2018

本研究通过模拟得到了4个湖泊所在区域的氮沉降通量值,并与已有文献[10, 15-16, 19]中的观测值进行对比,结果如表4所示。

表4 2018年青海湖、乌梁素海、东湖、太湖沉降通量模拟值与观测值比较Table 4 Comparison between simulated and observed nitrogen deposition fluxes of Lake Qinghaihu, Lake Ulansuhai, Lake Donghu, and Lake Taihu in 2018

尽管由于观测站点位置、观测时段、观测方法和数值模拟方法中沉降速率参数的计算等因素,模拟值与观测值存在一定差异,但沉降通量总体上具有可比性。其中,青海湖的湿沉降模拟值低于文献观测值,因Zhang等[15]的研究关注整个青海湖流域的氮沉降而将采样站点设在铁卜加草原改良实验站,该站点距离青海湖岸边约11 km,靠近西宁市区,周边排放源密集,空气中含氮污染物含量较高。本研究中湖泊沉降通量模拟值仅考虑青海湖面,由湖面所在区域网格计算平均值得到,而青海湖水面面积广阔,湖面上方含氮污染物浓度低于西宁市区,因此造成模拟值低于文献中的观测值。针对太湖不同年份和季节的大气氮沉降研究较多,例如,2002年7月—2003年6月湿沉降总氮为28.07 kgN/hm2[51],2015年1—12月干、湿沉降总氮分别为50.12、11.42 kgN/hm2[52],2014年3月—2016年2月年均干湿沉降总氮为64.8 kgN/hm2[11],即太湖年均干沉降总氮为8.21～50.12 kgN/hm2、湿沉降总氮为8.42～28.07 kgN/hm2,本研究模拟值也位于上述范围内。

4个典型湖泊大气氮沉降通量的季节模拟结果如图4所示。4个湖泊中,东湖与太湖均处于发达地区,人为排放源密度显著高于青海湖与乌梁素海。东湖面积较小,属于城市中心的城中湖,而太湖横跨两省五市,周边以森林与湿地公园为主,因此,东湖更易受到城市人为排放源的影响,氧化态氮与还原态氮沉降通量都明显高于太湖。夏冬两季东湖与太湖的总氮沉降通量较高,东湖分别为18.42 kgN/hm2与15.89 kgN/hm2,占全年29.6%和30.3%,太湖分别为6.06 kgN/hm2与6.21 kgN/hm2,占全年29.5%和27.0%。频繁的农业活动、高温天气导致电力需求增加、及氨挥发效率的增加是导致东湖和太湖夏季氮沉降通量偏高的主要因素。

图4 2018年4个典型湖泊还原态氮(REDN)与氧化态氮(OXN)沉降通量的季节模拟结果Fig.4 Simulated seasonal deposition fluxes of REDN and OXN in the four typical lakes in 2018

同时,2个湖泊冬季氮沉降通量与夏季相当,太湖冬季氮沉降通量甚至高于夏季,可能是因冬季大气臭氧浓度相对较低(见图2),大气呈现出较低的氧化性,导致湖区大气中还原态氮的累积,造成湖泊冬季氮沉降通量偏高。东湖夏季和冬季还原态氮沉降通量较高,分别为10.11、11.98 kgN/hm2,夏季氧化态氮沉降通量较高,为8.31 kgN/hm2。太湖夏季和冬季还原态氮沉降通量均较高,分别为2.49、2.34 kgN/hm2,氧化态氮沉降通量也较高,分别为3.57、3.87 kgN/hm2。另外,太湖临近长三角城市群,冬季能源消耗量相对较高,氮氧化物排放量增加,也是太湖冬季氧化态氮沉降通量较高的原因之一。4个湖泊还原态氮部分分别占57%、55%、60%、38%,除太湖每个季节还原态氮沉降通量都高于氧化态氮沉降通量外,其余湖泊都是还原态氮沉降通量占主要部分。在中国大部分地区,由于农业区广泛存在,总氮沉降通量通常由还原态氮主导[53-54],而太湖由于临近长三角城市群,工业源和交通源的氮氧化物排放较为集中,因此该湖区氧化态氮沉降通量比例较高。

4个湖泊中,青海湖与乌梁素海位于西北和北部偏远地区,全年的大气氮沉降通量较低。青海湖夏、秋两季氮沉降总通量分别为0.75、1.15 kgN/hm2,高于春、冬两季(分别为0.42、0.57 kgN/hm2),乌梁素海冬季氮沉降通量较低,仅为0.66 kgN/hm2,夏季较高,为1.18 kgN/hm2。青海湖与乌梁素海夏季沉降通量分别占全年的26.0%和28.6%,而冬季氮沉降通量分别占全年的19.7%和16.0%,除了氮氧化物排放源较少外,也与我国西北和北部地区冬季降雨较少有关。4个湖泊全年的氮沉降通量中,以湿沉降通量为主,分别占86.4%、56.3%、62.0%、73.2%;除了乌梁素海冬季干沉降通量略高于湿沉降通量外,其余湖泊在4个季节也都以湿沉降通量为主。

3.2 大气氮沉降对湖泊污染的影响分析

根据水利部太湖流域管理局发布的2018年太湖水健康状况公报[31]中公布的数据,2018年全年太湖流域入湖河流输入的氮总量为3.96×107kg。基于本研究中CMAQ模式模拟结果,太湖2018年总氮沉降通量为21.67 kgN/hm2,结合太湖面积计算可得2018年太湖大气沉降来源的氮负荷约为5.067 ×106kg,约占总氮负荷的11.3%。在此估算中,来源于入湖河流与大气沉降之外的氮负荷忽略不计。同理,可计算得到青海湖、乌梁素海与东湖2018年总氮沉降通量分别为3.02、4.38、64.25 kgN/hm2,乌梁素海与东湖的大气沉降来源的氮负荷比例分别为0.9%、22.3%(如表5所示)。因青海湖全年总氮负荷数据未见报道,氮负荷比例暂无法计算。

表5 4个典型湖泊及中国其他湖泊的大气沉降氮负荷Table 5 Nitrogen load from atmospheric depositions of four typical lakes and some other lakes in China

由表5可知,青海湖和乌梁素海的年大气氮沉降总量远低于东湖和太湖,乌梁素海大气氮沉降对总氮负荷贡献比例为0.9%,可见大气氮沉降对污染源排放密度低的地区的湖泊不会造成显著影响。根据其他研究湖泊氮沉降的文献报道,大气氮沉降对于湖泊氮负荷的贡献比例约为6.14%～48.8%[4,11,13,55-57],与本研究结果具有可比性。表5中,各研究结果呈现出一定程度的差异,主要原因为:观测研究多采用在湖泊岸边或湖心岛设置采样器收集降尘、雨水样品并对氮沉降量进行分析测算,采样点的设置以及样品采集过程本身的不确定性会导致误差;本研究采用数值模拟的方式测算湖泊的氮沉降,模式本身以及排放清单的不确定性也会使结果产生偏差。

3.3 大气氮沉降模拟的敏感性分析

图5 沉降速率系数与排放速率系数调整前后典型湖泊氮沉降通量模拟值对比Fig.5 Comparison of nitrogen deposition fluxes of typical lakes before and after adjusting the deposition and emission rate coefficients

氮元素沉降速率降低20%使得4个湖泊全年总氮沉降通量降低3%～11%,氧化态氮与还原态氮降低比率与总氮沉降通量接近一致;源排放速率降低20%可使得4个湖泊全年氮总沉降通量下降11%～22%,其中6月份氧化态氮沉降通量4个湖泊平均下降17.3%,还原态氮沉降通量4个湖泊平均下降18.7%,12月份氧化态氮沉降通量4个湖泊平均下降12.2%,还原态氮沉降通量4个湖泊平均下降22.4%。可见,氮沉降通量对于源排放速率的敏感性要高于对于沉降速率的敏感性,且源排放速率降低对于还原态氮沉降通量的影响要高于氧化态氮,具体原因需要未来进一步的探究。从大气沉降的角度揭示了氮排放源的削减对于湖泊氮污染防治的重要意义,同时也表明,CMAQ沉降通量的模拟较大程度依赖于含氮污染物的排放清单的准确性。

4 结论与展望

本研究对2018年中国4个典型湖泊的大气氮沉降通量进行了数值模拟,与现有观测值进行对比并分析沉降特征,估算大气氮沉降对湖泊污染的影响,通过敏感性分析识别影响氮沉降模拟过程的关键影响因子。所得结论如下:

(1)青海湖、乌梁素海、东湖与太湖全年氮沉降通量分别为3.02、4.38、64.25、21.67 kgN/hm2,大气氮沉降贡献的氮负荷分别占乌梁素海、东湖与太湖总氮负荷的0.9%、22.3%、11.3%,东湖与太湖受大气氮沉降的影响程度要显著高于青海湖和乌梁素海。

(2)4个湖泊的氮沉降通量有较显著的季节特征,夏季氮沉降通量都高于全年平均水平,分别占到全年的26.0%、28.6%、29.6%、29.5%。东湖冬季还原态氮沉降通量较高,为11.98 kgN/hm2;夏季氧化态氮沉降通量较高,为8.31 kgN/hm2。4个湖泊全年的氮沉降通量中,以湿沉降通量为主,分别占86.4%、56.3%、62.0%、73.2%。敏感性分析的实验结果表面,含氮污染物的沉降速率对湖泊氮沉降通量模拟的影响不显著,而含氮污染物的源排放速率的削减会使得氮沉降通量明显降低,源排放速率降低20%,4个湖泊总氮沉降通量下降11%～22%。因此,可通过完善含氮污染物排放清单等方式进一步提升CMAQ模式对于氮沉降通量模拟的准确性。同时,发达地区湖泊的大气氮沉降污染问题值得关注,且需要通过削减氮排放源等手段来进一步分析减少大气氮沉降对降低湖泊总氮负荷的作用。

本研究仅使用WRF-CMAQ模式对2018年我国不同地区的4个典型湖泊进行模拟,覆盖范围较广,选取了相对普适的模拟参数方案,未来还可针对不同年份开展对比研究,并针对某一湖泊,选取特定的模拟参数方案,进行精度更高的模拟研究,从不同的时间和空间尺度上深入探索湖泊的大气氮沉降特征。