基于欧拉-拉格朗日方法的气载核污染物长距离传输沉降模拟研究

张 彦 乔清党 杨端节 郭 猜 刘金凯

（生态环境部核与辐射安全中心 北京 100082）

核爆炸或重大核泄漏事故产生的核污染物由于可能对人体健康和环境带来重大危害和影响，其在大气中的长距离传输和沉降过程，一直是各国核应急工作者关注的重点。通过数值计算方法，开展模拟研究是预测和重现核污染物在大气中的传输路径及对污染后果进行定量评估的重要手段。

目前，在中长距离核污染物扩散模拟和预报的数值方法上，应用比较成熟和广泛的主要是拉格朗日方法和欧拉方法两类，分别基于湍流统计理论和梯度输送理论。各国核安全监管当局和科研机构，在采用拉格朗日或欧拉方法的基础上陆续开发了多种数值模拟软件系统，如欧盟的RODOS（Real-time On-line Decision Support）［1-2］、日 本 的SPEEDI/WSPEEDI（System for Prediction of Environmental Emergency Dose Information/Worldwide version of System for Prediction of Environmental Emergency Dose Information）［3］、美 国 的ARAC/NARAC（Atmospheric Release Advisory Capability/National Atmosphere Release Advisory Center）［4］和HYSPLIT（Hybrid Single-Particle Lagrangian Integrated Trajectory Model）［5］、我国的RADCON（Radioactive Consequence Assessment System for Overseas Nuclear Accident）［6］、PARMODEL（Particle diffusion Model）［7］等。两种方法都有各自的优缺点，拉格朗日方法更适用于对统计平稳和均匀湍流条件下的扩散问题描述，但对各种大气物理过程及干、湿沉降考虑都比较简单；欧拉方法能够模拟更复杂的物理过程，但是容易出现网格点上虚假扩散等问题［8］。欧拉-拉格朗日方法则通过两种数值模拟方法的结合，能够进一步提高模型的适用性和应用范围，国外对此研究相对较多。Nikolai Talerko应用基于LEDI（Lagrangian-Eulerian DIffusion model）大气传输模型，重构了乌克兰地区由于切尔诺贝利核事故产生的放射性污染随时间变化，给出了131I的放射性污染图形［9-10］；Stein和Cohen等［5］将基于拉格朗日方法（HYSPLIT）和欧拉方法的全球模型（Global Eulerian Model，GEM）相嵌套，开展了欧拉-拉格朗日方法模拟研究，得到了较好的模拟结果。近年来，欧拉-拉格朗日方法在水体及密闭空间内污染扩散等方面也得到了一定的应用研究［11-12］。

本文基于JRODOS（Java Real-time On-line Decision Support）系 统 MATCH（Multi-scale Atmospheric Transport and Chemistry）模块的欧拉-拉格朗日方法，针对我国东部沿海某核电厂假想的核事故泄漏，开展长距离气载核污染物传输扩散数值模拟研究，进一步探讨欧拉-拉格朗日方法在国内核事故污染物长距离输运模拟中的适用性。

1 数值模拟软件系统简介

欧盟的RODOS软件系统是在欧洲得到广泛应用的综合性实时在线场外核应急决策支持系统，它可以在整个事故阶段和各个距离范围对事故后果、防护行动和对策作出评价、分析和预测［1-2］。近年来，RODOS系统广泛使用的版本为JRODOS，是Java版本的RODOS系统，可以实现Windows和Linux版本的统一运行，用户界面和体验性得到极大增强。其最初版本于2009年12月在德国KIT（Karlsruhe Institute of Technology）运行中心完成测试，2010年4月发布第二个运行版本。2011年3月，针对日本福岛核电厂定制了本地化的JRODOS，并用于开展福岛核事故的核污染物大气扩散和沉积模拟计算工作［13］。2019年生态环境部通过中欧合作项目，与欧盟达成相关协议，获得了JRODOS最新版本的使用权，并将其作为生态环境部（国家核安全局）核与辐射事故后果评价工作的技术支撑工具之一。

在长距离核污染扩散方面，欧盟将瑞典气象水文研究所（Swedish Meteorological and Hydrological Institute，SMHI）开发的多尺度大气输送和化学模型MATCH（Multi-scale Atmospheric Transport and Chemistry）［14］嵌入RODOS系统中用于相关模拟工作。MATCH模型在最初开发时完全采用了欧拉方法，主要应用于空气质量预报、化学污染物传输扩散等多个领域。1998年，Langner等［15］将其应用在放射性物质长距离扩散研究中，通过在放射性源项处理上使用拉格朗日方法进行初始化，将其扩展为基于欧拉-拉格朗日方法的数值模型。

欧拉-拉格朗日方法模拟的基本思路是：在污染物释放的初始阶段，采用拉格朗日方法来描述其扩散传输过程，当传输经过一定距离且污染物扩展到足够多的计算网格之后，采用欧拉方法对后续扩散进行模拟。对此，MATCH模型在源项初始扩散10 h内采用拉格朗日方法，10 h后切换为欧拉方法，同时将源项转换为欧拉方法需要的浓度场。模型采用的拉格朗日方法仅用来处理平流和水平扩散，在垂直扩散和干、湿沉降方法上仍按照欧拉方法来处理［15］。MATCH模型针对放射性污染物粒径进行了分型，赋予其不同的重力沉降速度，能够用来模拟核爆炸烟云扩散过程［16］，但在JRODOS系统 引入的MATCH模块中，不提供核爆后果评价计算功能，在核电厂核事故后果评价计算中，将放射性污染物看作气溶胶粒子，按照相关经典干、湿沉降理论处理［15］。2014年，Kovalets等［17］应用JRODOS系统的MATCH模块对日本福岛核电厂核事故污染扩散进行了相关研究，取到了较好的结果。

2 核事故模拟输入设置

2.1 事故源项设计

2011年3月11日，日本福岛核电厂因地震海啸导致发生7级核事故，大量核污染物泄漏到大气中，事后日本及国际上各方都对泄漏数据进行了推算及反演研究。根据联合国原子辐射影响科学委员会（United Nations Scientific Committee on the Effects of Atomic Radiation，UNSCEAR）在2013年报告书中的权威评估，福岛核事故的泄漏量：131I为100~500 PBq，137Cs为6~20 PBq，约为切尔诺贝利核事故泄漏量的1/10和1/5；日本先后撤离居民约88 000人，成人在撤离前及撤离中所受有效剂量平均不到10 mSv，从数据上来说对全球人和环境造成的辐射危害相对有限，但对当地居民生产和生活带来了不可逆的后果，其影响仍然超出了社会和公众的接受限度［18］。

本次模拟的源项设计参考福岛核事故的最大泄漏量，假想我国东部沿海地区某核电厂，在出现某种不可抗力的意外情形下，发生超设计基准严重事故，导致核污染物释放到大气中，采用JRODOS系统的MATCH模块对其长距离传输、扩散、沉降过程开展数值模拟研究。

JRODOS系统的MATCH模块中，可对10余种核素进行模拟计算，包括：88Kr、133Xe、135Xe、131I、132I、133I、135I、88Rb、89Sr、90Sr、95Zr、132Te、140Ba、134Cs、137Cs等核素，对不同核素按照其各自半衰期分别计算，为简化计算未考虑各核素的衰变链。结合福岛核事故核污染物释放情况，本次模拟计算选取131I、137Cs两种典型核素作为释放源项，并假设131I释放总量为500 PBq，即5×1017Bq，137Cs释放总量为20 PBq，即2×1016Bq，131I释放总量为137Cs的25倍。

2.2 模型参数初始化设置

MATCH计算需要大尺度数值天气预报数据作为输入气象场，支持NCEP（National Centers for Environmental Prediction）、ECMWF（European Centre for Medium-Range Weather Forecasts）等多种国际通用气象数据格式，通过JRODOS内部工具进行数据转换后，可生成计算所需的气象场，传输扩散计算的分辨率由气象场的空间分辨率同步确定。本次模拟计算随机选取了2021年1月25日18时（UTC时间）到2021年1月29日21时（UTC时间）共99 h的NCEP全球GFS预报资料作为背景气象场［19］，其中初始时刻25日18时数据为再分析数据，其后每隔3 h产生一个预报数据，逐次递增至99 h。数据内容覆盖各等压面（从地面到高空0.4 hPa高度），主要包含位势高度、温度、位温、风、海平面温度、降水和对流顶层气压等物理量。选取的气象数据空间分辨率为0.5°，模拟区域以我国东部某核电厂（30.5°N，120.9°E）为中心，网格范围设置为100×100，模拟范围为95°E~145°E，6°N~56°N，覆盖我国东部、东南部及东亚、东南亚部分地区，水平扩散网格与气象场网格保持一致，垂直方向通过线性插值处理，计算结果的输出显示，需要与JRODOS系统的GIS（Geographic Information System）地图匹配，本次计算设置为墨卡托投影，计算结果按照墨卡托投影进行变换后输出。

核事故源项释放时间设定在2021年1月25日19时（UTC时间），假设释放高度为1 km，释放时间总计为3 h。系统需要输入源项的平均释放速率，根据假定释放总量推算，设置释放平均速率131I约为4.6×1013Bq·s-1，137Cs约为1.85×1012Bq·s-1。

3 模拟结果与分析

3.1 长距离传输沉降总量模拟结果分析

由于源项设置中131I释放总量为137Cs的25倍，模拟计算结果表明，在地面沉降份额中131I所占比例达到96%左右，137Cs为3%左右。图1为本次核事故模拟中131I在长距离传输中的地面累积沉降量，包括干沉降和湿沉降总和，选取的4个时刻分别为26日00时UTC、26日09时UTC、27日00时UTC、29日21时UTC，即释放结束后2 h、11 h、26 h和95 h的地面累积沉降量分布图。

分析图1结果可以看出，图1（a）、（b）中的地面沉降分布明显具有拉格朗日方法计算的结果特征，即与空中单一污染烟团相对应，呈现出规则的矩形叠加分布，同时也体现出早期污染烟团传输距离较短，水平扩散范围较小的明显特征；其后，从图1（c）、（d）可以看到，随着风向的变化以及模型转为欧拉计算方法，地面沉降图形逐步表现为大范围的蘑菇状扩散，沉降图形由从事故点向东南方向扩展，转为从东南向西南方向延伸，且在向西南主方向延伸的同时，保持着在各个方向上的随机扩散形态，明显具有欧拉方法中各网格点独立扩散的特征；最终在计算终止时，地面沉降图形从我国华东地区延伸至东海、台湾海峡呈带状分布，覆盖整个华南沿海地区，包括台湾和海南岛全部，并一直延伸至越南北部，影响面积东西方向达到3 000 km左右，南北方向达到2 000 km左右。

图1 某核电厂核事故131I长距离传输地面沉降模拟结果（Bq·m-2）(a) 1月26日00时UTC，(b) 1月26日09时UTC，(c) 1月27日00时UTC，(d) 1月29日21时UTCFig.1 Surface deposition simulation results of 131I long distance transport of nuclear accident in a nuclear power plant (Bq·m-2)(a) 00:00 UTC, January 26, (b) 09:00 UTC, January 26, (c) 00:00 UTC, January 27, (d) 21:00 UTC, January 29

在沉降量级上可以看到，最大沉降量在107~108Bq·m-2区间，出现在核事故发生早期，距离核电厂较近，在核电厂下风东南方向，影响范围约占4个网格面积，约10 000 km2；在随后的沉降图形上，由蘑菇状分布逐渐转为明显的带状分布，沉降量级逐次递减至10 Bq·m-2，范围逐渐扩大，在西南方向下游出现103~104Bq·m-2和102~103Bq·m-2的局部中心点。由此说明，在核事故早期，应重点关注距事故地点下游方向的近距离污染情况，中期（事故发生30 h后）应关注在距事故地点远距离外可能出现的沉降中心点。

图2为本次核事故模拟中137Cs在长距离传输中的地面累积沉降量，包括干、湿沉降总和。总体看，137Cs地面沉降图形与131I极为相似，但由于释放总量仅为131I的1/25，导致其在地面沉降贡献占比中仅为3%左右，在各个时刻的沉降图形覆盖面积也比131I明显偏小；137Cs早期最大沉降量级在106~107Bq·m-2区间，出现时间和地点与131I沉降图形非常接近，在后期下游的局部中心点，沉降量级均比131I约低一个量级，中心区域范围也明显缩小。模拟结果表明：131I在污染范围和量级上的贡献远超137Cs，与最初的释放源项份额差异有明显关联。

图2 某核电厂核事故137Cs长距离传输地面沉降模拟结果（Bq·m-2）(a) 1月26日09时UTC，(b) 1月29日21时UTCFig.2 Surface deposition simulation results of 137Cs long distance transport of nuclear accident in a nuclear power plant (Bq·m-2)(a) 09:00 UTC, January 26, (b) 21:00 UTC, January 29

3.2 湿沉降模拟结果分析

本次模拟计算采用的NCEP气象预报数据表明，在核污染烟云传输路径上有较长时间的大范围降水过程。图3为本次核事故131I地面湿沉降及降水预报模拟结果，137Cs与之相类似，为便于分析比较，选取截图时间与图1时间点相同；图4为131I在29日21时UTC地面湿沉降与干沉降模拟结果的比较图。

图3 某核电厂核事故131I地面湿沉降(a)（Bq·m-2）及降水预报(b)（kg·m-2）模拟结果(a1、b1) 1月26日00时UTC，(a2、b2) 1月26日09时UTC，(a3、b3) 1月27日00时UTC，(a4、b4) 1月29日21时UTCFig.3 Simulation results of 131I surface wet deposition (a) (Bq·m-2) and precipitation forecast (b) (kg·m-2) of nuclear accident in a nuclear power plant(a1, b1) 00:00 UTC, January 26, (a2, b2) 09:00 UTC, January 26, (a3, b3) 00:00 UTC, January 27, (a4, b4) 21:00 UTC, January 29

从图3（a）中在4个时间点上的截图可以看出，131I的地面湿沉降在范围和量级上与图1总沉降结果极为接近，图4也表明，干沉降分布范围相对湿沉降小很多，且最大值与湿沉降结果相差两个量级左右，说明在本次模拟的沉降结果中，湿沉降占主要地位；在图3（b）中对应时间的降水预报结果上也可以看到，在核污染物传输下游一直有持续性的降水出现，在图3（b4）中，台湾海峡和越南北部附近还出现了局部的较高降水中心区，与图3（a4）中的湿沉降局部中心点呈现出对应关系。以上分析可以看出，本次模拟中降水对核污染物的清除发挥了重要作用，特别在核事故早期降水的清除作用，也是导致核污染物远、近距离沉降相差4~5个量级的主要原因。

图4 某核电厂核事故131I地面湿沉降(a)（Bq·m-2）及干沉降(b)（Bq·m-2）1月29日21时UTC模拟结果比较Fig.4 Comparison of simulation results of 131I surface wet deposition (a) (Bq·m-2) and surface dry deposition (b) (Bq·m-2) of nuclear accident in a nuclear power plant at 21:00 UTC on January 29

3.3 实时气象场对比分析

本次模拟计算采用的是气象预报数据，与应急响应时开展核事故后果评价工作相一致。为进一步验证气象预报数据和模拟结果的准确性，选取了模拟计算时段的天气分析图进行对比验证。图5为我国中央气象台在模拟释放后前两天的部分实时地面天气叠加雷达降水分析图［20］。

图5 中央气象台实时地面天气叠加雷达降水分析图(a) 1月25日21时UTC，(b) 1月26日00时UTC，(c) 1月26日09时UTC，(d) 1月27日00时UTCFig.5 Real-time surface weather superimposed radar precipitation analysis chart of Chinese Central Meteorological Station(a) 21:00 UTC, January 25, (b) 00:00 UTC, January 26, (c) 09:00 UTC, January 26, (d) 00:00 UTC, January 27

从图5可以看出，25日21时UTC，即在核事故污染物释放过程中，事故地点以西北风为主，风向下游及我国东部沿海地区有明显的实时降水出现，后续几天内，东南沿海雷达图显示沿岸站点一直有小规模降水，与NCEP气象预报数据较为一致；在计算时间段内，我国华东、华南地区一直处于地面高压前部控制，等压线较为稀疏，风速较小，高压系统稳定且逐步向东南移动，东南沿海及台湾海峡以东北风为主，因此污染烟云在传输过程中，呈现出先向东南再转向西南方向缓慢移动扩散为主的路径，与地面天气实况相一致。

图6为我国中央气象台在模拟释放后前两天内的部分高空实时天气分析图［20］。本次模拟设置的源项释放高度为1 km，故选取了925 hPa（约800 m）和 850 hPa（约1 500 m）高空天气图开展对应分析。1月26日00时UTC，即释放结束2 h后，925 hPa等压面上，事故地点处于高压前部，受弱冷平流影响，主导风向为西北风；850 hPa等压面上，事故地点处于高低压之间鞍型场控制，弱冷平流影响比925 hPa更小，主导风向为偏西风。整体看，污染烟云处于两层之间，由于两层等高线都较为稀疏，风速较小，系统稳定少动，垂直运动较小，不利于污染烟团快速向远距离传输，因此，在降水作用下事故地点近距离容易形成沉降最大值。1月27日00时UTC，即释放结束26 h后，烟云主体已经扩散至我国东海、台湾海峡大部，该处气象场在925 hPa等压面上处于高压底部，受弱冷平流影响，主导风向为偏北风；850 hPa等压面上受南部高压整体控制，主导风向为东北风及东风，受弱暖平流影响。两层高空等高线仍然较为稀疏，风速较小，系统移动缓慢，不利于垂直扩散和长距离快速传输。因此，污染烟云在此时逐渐缓慢向我国南部移动，最终沿着台湾海峡形成从东北-西南走向的狭长带状分布沉降图，数值模拟结果与实况天气变化趋势较为吻合。

图6 中央气象台925 hPa (a)及850 hPa (b)实时高空天气图(a1、b1) 1月26日00时UTC，(a2、b2) 1月27日00时UTCFig.6 Real-time upper-air chart of 925 hPa (a) and 850 hPa (b) of Chinese Central Meteorological Station(a1, b1) 00:00 UTC, January 26, (a2, b2) 00:00 UTC, January 27

3.4 地面剂量率场模拟结果分析

为进一步分析核污染沉降对人群健康和环境的影响，模拟结果还给出了地面总剂量率场（包括131I和137Cs总贡献）分布图（图7）。由于131I释放总量为137Cs的25倍，模拟结果显示，131I对地面剂量率的贡献大于99%，与总剂量率场分布结果基本完全一致，表明本次模拟中131I对人和环境的影响比137Cs更加重要。

图7的模拟结果表明，地面剂量率的最大量级在10-5~10-4mSv·h-1区间，即为10~100 nSv·h-1区间，和天然本底水平相当，在核事故源项释放结束11 h后出现在我国东部海域；其后最大值下降为10-7~10-6mSv·h-1区间，即0.1~1 nSv·h-1区间，远小于天然本底水平，28日00时之后地面剂量率场量级进一步下降，低于最低绘图区间值。总体看，在各时间点上，地面剂量率的出现范围也小于总沉降图形范围，且主要影响范围为海上无人区，量级上接近或低于天然本底水平，表明本次模拟核事故的污染沉降对人体健康和环境影响较小。

图7 某核电厂核事故长距离传输地面剂量率场（mSv·h-1）模拟结果(a) 1月26日09时UTC，(b) 1月26日19时UTC，(c) 1月27日00时UTC，(d) 1月28日00时UTCFig.7 Surface dose rate field (mSv·h-1) simulation results of long distance transport of nuclear accident in a nuclear power plant(a) 09:00 UTC, January 26, (b) 19:00 UTC, January 26, (c) 00:00 UTC, January 27, (d) 00:00 UTC, January 28

4 结语

通过本次数值模拟及结果分析，可以得到以下初步结果：

1） JRODOS系统采用的欧拉-拉格朗日模型方法能够给出核污染物长距离传输和沉降的主要特征，本次模拟99 h污染扩散耗时约12 min，计算效率较高，结果可视化程度较好，模拟结果与实况天气图总体趋势上较为符合，可以为我国核事故后果评价及应急决策提供辅助参考。

2） 在本次模拟的气象条件下，在发生类似福岛核电厂泄漏的7级核事故时，其对我国东部、东南沿海地区及台湾、海南岛等地的影响，地面剂量率场最大值在10~100 nSv·h-1区间，后续下降为0.1~1 nSv·h-1，相当或低于天然本底水平。

3） 本次模拟结果表明，在天气形势不利于污染物向长距离快速传输的情况下，降水导致的湿沉降将对污染物的快速清除起到决定性作用，同时传输路径上的局部强降水中心将导致地面沉降出现局部中心高值，应引起重点关注。

作者贡献声明张彦：负责文章编写，开展模型实验并分析结果；乔清党：参与结果分析讨论，负责文章校核；杨端节：提出修改意见、审核论文；郭猜：参与文献调研及结果分析讨论；刘金凯：参与结果分析讨论，协调办理相关发表流程。