气象扩散耦合模型对经典示踪试验ETEX-1的模拟分析

郭 欢,张俊芳,李鸣野,黄 莎,赵多新

(中国辐射防护研究院 核环境模拟与评价技术重点实验室,太原 030006)

核事故发生后,气载放射性污染物在大气中迅速扩散,短时间内可能造成大范围的环境影响,若能精确预测污染物的扩散规律,就可为及时制定应急响应决策提供参考[1]。在评估污染物在大气中的扩散过程时,单纯依靠观测资料存在数据有限、无法预报未来扩散趋势等短板问题,而数值模拟因其不受试验条件和天气条件限制等优点成为开展大气扩散研究的一种常用手段[2]。大气扩散主要受大气边界层中的湍流控制,而湍流存在随机性,通常采用统计学方法将野外示踪试验数据与数值模式的模拟结果作比对来描述其对扩散的影响[3-4]。陈晓秋[5]利用大亚湾核电厂的示踪试验数据对随机游走模式的数值模拟能力进行了验证;胡二邦等人[6]基于福建惠安核电厂址的野外示踪试验数据对高斯烟羽模式在复杂滨海地形下的模拟进行了分析研究;王博[7]以野外示踪试验数据为参考,对CFD的大气扩散特征模拟能力进行了验证;朱好等人[8]利用湖南某厂址的野外示踪试验数据对拉格朗日烟团模型CALPUFF在复杂地形下的近场适用性开展了研究。

本文采用WRF V3.5中尺度气象模式和FLEXPART-WRF V3.1粒子扩散模式,并根据WRF版本,调节WRF输出数据中的默认气象要素以及调节模式输出时空分辨率、模拟区域、物理参数化方案等参数信息,改动接口程序以适应V3.1版本粒子扩散模式的输入需求,利用国外经典示踪试验ETEX-1数据集和由丹麦国家环境研究所(NERI)开发的MVK工具中的统计方法,对试验期间的情景进行了数值再现,并将模拟结果与相关的气象观测数据及现场浓度数据作比对,对耦合模式模拟效果进行了验证分析。

1 模式设置与研究方法

1.1 ETEX试验概况

欧洲长距离示踪试验(ETEX)[9]包括在1994年10月和11月期间开展的两次试验,根据气象背景统计分析,中欧上空出现西风气流时大多数欧洲国家会受到示踪羽流的影响,天气条件最有利于试验开展,因此释放点定于法国西部的蒙特菲尔。为了避免交叉污染和样本污染,第一次释放PMCH(全氟甲基环己烷)作为示踪物,第二次释放PMCP(全氟甲基环戊烷),且两次试验的时间间隔超过2周。试验中烟羽迁移范围达2 000余km,在释放后72小时内每间隔3小时进行采样,使用遍布欧洲17个国家的168个地面监测站点的采样网收集空气样本,两次实验采集到约9 000份样本。示踪数据库中汇集了释放特征信息、地面和高空的常规与附加气象观测、地面和高空浓度等数据。第一次释放(ETEX-1)期间,天气状态平稳,没有出现高压或低压中心,也没有延伸的脊或槽靠近释放点,无锋面系统经过释放地点,在稳定的强西风和西南风气流引导下,示踪羽流被输送到整个欧洲。第二次试验(ETEX-2)期间受释放初期锋面系统过境的影响,天气条件恶劣,烟羽并未按照预期覆盖所布监测点,对试验的开展和空气样本数据质量都产生了极大不利影响,国际上通常认为第一次试验较为成功,因此本文选择第一次试验(ETEX-1)期间的情形作模拟分析。

图1为示踪试验期间的站点分布,包括欧洲168站和法国西北部81站,其中,欧洲168站既进行气象观测又收集浓度样本,法国81站仅进行常规气象观测。

红:释放点;蓝:欧洲168站;绿:法国81站。

1.2 WRF 模式参数设置

根据示踪试验时段,选用欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts,ECMWF)的再分析数据作为WRF模式的初始场和边界条件,对试验期间的72 h进行模拟,即1994年10月23日12时至1994年10月26日12时(UTC),数据空间分辨率为1°×1°,地形数据采用美国地理学会(USGS)的10′分辨率地形数据。对模拟的欧洲地区进行网格划分,模拟中心位置设为49°N,12°E,采用两层嵌套,第一层(d01)网格分辨设定为30 km,模拟区网格数为111×111,第二层(d02)网格分辨设定为10 km,模拟区网格数为231×231。在对物理参数化方案的设置上,微物理参数化方案选择WSM3,长波和短波辐射方案分别选择RRTM和Dudhia,边界层参数化方案选择MYNN2.5,陆面过程方案选NOAH。

1.3 FLEXPART-WRF模式参数设置

FLEXPART-WRF由驱动场气象模型WRF和主体扩散模型FLEXPART组成。WRF模型以原始全球数值预报产品作为输入,经过模式动力降尺度后可获取更高时空分辨率的高质量气象场,从而为扩散模型提供更精确的气象初始场。FLEXPART模式由挪威大气研究所开发,是一种典型的拉格朗日粒子扩散模式,通过设置释放大量随机粒子,跟踪其在湍流中的运动来获取粒子位置和扩散浓度[10]。FLEXPART在模拟过程中引入边界层参数化方案来反映对流和湍流运动特征,因此可以更准确地描述粒子在大气中的运动过程。两种模式经过耦合,通过改善气象场模拟效果可提高扩散模拟精度。但耦合模型也存在一定的局限性,其描述的是几百至几千km的大范围扩散过程,对于几km至几十km范围的扩散模拟,还需对模型做更深入的研究和改进。

WRF输出为FLEXPART-WRF提供必要的输入,WRF输出具体变量见表1。

表1 FLEXPART-WRF输入所需的WRF输出变量

扩散模拟以示踪试验条件为参考,FLEXPART-WRF模式中,释放点为法国蒙泰尔菲(48°03′30″N,2°00′30″W),示踪物为PMCH,总释放量340 kg,释放高度为离地8 m,起始释放时间设为1994年10月23日16时(UTC时间),连续释放12小时。

1.4 验证方法

NERI以Hanna等[11]提出的定量统计参数为基础,研发了一套模型验证比对工具MVK(Model Validation Kit),其中的定量分析程序(Boot)包含全面的ASTM(American Society for Testing and Materials)模型验证评价方法,通过将试验观测值和模型预测值进行定量统计分析,可以有效评价模型的准确性。统计指标包括平均值(Mean)、平均偏差(Bias)、标准偏差(Sigma)、偏差分数(FB)、几何平均偏差(MG)、几何平均方差(VG)、标准化均方误差(NMSE)、相关系数(R)等。另外,美国国家海洋和大气管理局(NOAA)组织的DATEM项目专门用于大尺度大气扩散模式的有效性评估,评估中使用相关系数(R)、偏差分数(FB)、空间吻合度系数(FMS)、KS检验(Kolmogorov-Smirnov,KS)组合为一个模式评分指标RANK对模式有效性进行评价。本文结合MVK工具和NOAA对模式的评分指标评价方法,对耦合模式进行了ETEX-1试验的再现模拟分析,并将分析结果与美国空气资源实验室(ARL)官网公布的HYSPLIT模式对ETEX-1模拟结果的RANK指标作对比。相关计算公式如表2所示。

表2 指标说明

2 结果分析

2.1 风场模拟结果分析

在核事故应急预警中,浓度场是重要的评估对象,而风场作为大气扩散的驱动因子会直接影响浓度分布,因此首先对风场进行分析。图2表示不同时刻模拟的地面风场。ETEX-1试验覆盖面广,范围涉及欧洲多个国家和地区。与示踪试验相适应,模拟区域覆盖整个西欧以及部分东欧地区,且包含北大西洋的一部分、北海、地中海、波罗的海等海域,试验期间模拟区西部有来自北大西洋的海风,在地中海、北海附近多出现环流。在试验过程中的大多数时刻,释放点均处于北大西洋吹向陆面的海风之下,且受地形和背景天气形势的共同影响,气流在法国东部开始分流,一部分融入北海环流中,一部分流向地中海。定性分析结果认为风场较好地模拟了地面风场受海陆影响而产生的扰动变化情况。

图2 不同时刻的模拟地面风场

对欧洲168个常规地面气象站、法国西北部81个常规地面气象站的观测结果与模拟的地面风场分别进行对比,对比结果见图3～图4,法国站U、V风矢量的模拟值与观测值之间的相关性比欧洲站高,均方根误差比欧洲站低,除个别时刻外,相关系数和均方根误差随时间的波动幅度较小。统计分析全部站点数据后,得出U、V风分量的相关系数均值为0.43,均方根误差均值为3.4 m/s,总体上,风场模拟结果能够代表实际风场情况。由于扩散与风场直接相关,因此风场从一定程度上为扩散模拟奠定了良好的气象条件基础。

图3 U、V风分量相关系数随时间的变化

图4 U、V风分量均方根误差随时间的变化

2.2 扩散场模拟结果分析

扩散与风场形势密切相关,示踪剂被释放后,随着模拟时间的推移,示踪粒子跟随风场向远距离迁移扩散,且示踪物分布范围逐渐变广。图5(a)给出示踪粒子扩散所造成的近地面浓度分布情况,结合模拟的地面风场分布进行分析,从释放点到远距离的扩散过程中,粒子呈现主要往东北和部分向东南迁移的趋势,近地面过境浓度分布与模拟风场结果匹配较好,从烟羽覆盖区和采样点的分布上也可以看出,模式对烟羽的空间位置捕捉效果较好。图5(b)为示踪粒子在整个试验阶段的最大过境浓度分布,其向东北的走势与风场形势吻合度很高,且烟云覆盖区域与欧洲相应的最大地面浓度监测站点的分布近乎一致。此外,根据模式计算结果,模拟最大浓度出现在释放点附近,与实际相符,说明耦合模式模拟结果能较好地反映风场对气载污染物的聚集和扩散的影响。

图5 近地面浓度分布

统计浓度模拟值与168站浓度监测值之间的平均绝对偏差,以及各自浓度不为0的站点个数,如表3所示。其中,时长表示从释放开始每隔3 h排序,这是为了与3 h时间分辨的试验采样数据保持一致。结果表明:在前24 h,观测与模拟的平均绝对偏差呈现先增大后减小的趋势,最高偏差为1.37 ng/m3,出现在释放后的第15 h。第27 h至第36 h出现短暂偏差增大的过程,随后逐渐较小,趋于平稳。对于浓度不为0的统计结果,观测和模拟个数均逐渐增大,这是由于示踪粒子一直在运动迁移,必然逐渐覆盖新的区域。结合模拟浓度的空间分布情况,释放63 h后,烟羽主要部分超出了采样网捕捉范围,因此,中欧虽仍存在低浓度测量值,但模拟浓度不为0的个数相比观测明显减少。

表3 浓度平均绝对偏差和浓度不为0的点位数据统计

散点图和分位数图在一定程度上能反映模型预测情况。在散点图中,观测值和预测值成对出现,可以直观显示模型预测过高或过低的程度。在统计分析时,剔除观测数据中标记为“-0.88”和“-0.99”的数据,其中,“-0.88”表示该站点位置能够检测到浓度,但未能给出有效数据,“-0.99”表示该站未检测到浓度值。图6中,数据点多落于左上角区域,说明模式模拟值比观测值整体偏高。分位数图(图7)将每个观测和预测数据从最低到最高进行排序,可以反映模式测预测趋势,从图中也可看出,模式对浓度的模拟总体偏高。

图6 模型散点图

图7 模型分位数图

选取2个指标分析其随时间的变化,从时间分布的角度考察模式性能(如图8所示)。Bias表征预测均值和观测均值的差值,从图8来看,除个别时刻Bias为负值,在模拟的绝大多数时刻均为正值,到后期二者差距逐渐接近于0,说明模式对浓度的模拟整体偏高。NMSE在释放开始后第3 h至第18 h较大,这可能是由于在短时间内释放粒子数较多,模拟扩散不充分造成局部浓度积聚,从而导致观测与模拟间产生了较大差距,之后减小逐渐稳定趋于0。

图8 模拟浓度与观测值的Bias、NMSE随时间的分布

计算统计学指标汇总如表4所示。FB和MG均可反映系统误差,FB为无量纲数,常用于对多量级浓度的比较,个别观测或预测浓度出现极大值会对其产生较大影响,因而在观测和预测完全相异时FB仍有可能为0,而MG因采用对数法对高值不敏感,但观测或模拟值近0会对其产生显著影响。VG反映由系统误差导致的几何方差和非系统误差导致的几何方差之和,对于预测浓度和观测浓度多量级变化的情况,更宜评估MG和VG。NMSE为偏差统计参量,对观测和预测间的偏差表现敏感。理想状态下,数值模拟的MG、VG值为1,FB和NMSE为0,但由于自然环境中大气运动存在随机性,预测结果很难达到此标准。FMS表示浓度观测区域和预测区域重叠的百分比,反映了模式的空间模拟效果[12]。KS(Kolomogorov-Smirnov)定义为两个累积分布之间的最大差值,是唯一根据未配对数据计算的统计参量,可衡量模型对实际观测浓度分布的再现程度,任何两个分布之间的最大差异不能超过100%。Chang等[13]综合考虑试验情景复杂程度以及气象条件的随机波动等影响因素,在不要求模拟值和实测值在时间和空间一一对应的条件下,推荐了指标的可接受范围,即-0.3

3 结论与讨论

本文将中尺度气象模型WRF和拉格朗日粒子扩散模型FLEXPART-WRF适应性耦合后进行了ETEX-1试验的再现数值模拟,以试验数据集为参考,结合MVK统计学指标,将模拟结果与气象观测和浓度监测数据分别作对比分析。从扩散整体情况看,粒子扩散趋势与流场相匹配,近地面过境浓度覆盖大部分监测站点;模拟与观测之间的偏差、标准化均方误差等指标值基本符合国际标准,说明模式有效地模拟了试验期间的天气过程和区域流场及扩散特征。本次模拟得到RANK值为1.89,比ARL实验室官网HYSPLIT对ETEX-1模拟结果的RANK小0.07,二者非常接近,认为气象扩散模型相耦合可有效反映放射性物质的迁移扩散规律,可用于辅助评价核事故下气载放射性污染物扩散对环境的影响。

需要说明的是,为了适应对ETEX-1试验的再现需求,本文采用再分析气象数据作为初始数据,在实际应急时,模式也可利用实时预报数据开展预报分析。此外,在未来的研究工作中,还需重点关注WRF模式中物理参数化过程的优化分析,以及采用多种气象观测资料进行数据同化以期获取更精确的气象场,从而进一步改善扩散模式的模拟效果。