張 彥 喬清黨 楊端節 郭 猜 劉金凱
（生態環境部核與輻射安全中心 北京 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）模塊的歐拉-拉格朗日方法，針對我國東部沿海某核電廠假想的核事故泄漏，開展長距離氣載核污染物傳輸擴散數值模擬研究，進一步探討歐拉-拉格朗日方法在國內核事故污染物長距離輸運模擬中的適用性。
歐盟的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.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.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
通過本次數值模擬及結果分析，可以得到以下初步結果：
1） JRODOS系統采用的歐拉-拉格朗日模型方法能夠給出核污染物長距離傳輸和沉降的主要特征，本次模擬99 h污染擴散耗時約12 min，計算效率較高，結果可視化程度較好，模擬結果與實況天氣圖總體趨勢上較為符合，可以為我國核事故后果評價及應急決策提供輔助參考。
2） 在本次模擬的氣象條件下，在發生類似福島核電廠泄漏的7級核事故時，其對我國東部、東南沿海地區及臺灣、海南島等地的影響，地面劑量率場最大值在10~100 nSv·h-1區間，后續下降為0.1~1 nSv·h-1，相當或低于天然本底水平。
3） 本次模擬結果表明，在天氣形勢不利于污染物向長距離快速傳輸的情況下，降水導致的濕沉降將對污染物的快速清除起到決定性作用，同時傳輸路徑上的局部強降水中心將導致地面沉降出現局部中心高值，應引起重點關注。
作者貢獻聲明張彥：負責文章編寫，開展模型實驗并分析結果；
喬清黨：參與結果分析討論，負責文章校核；
楊端節：提出修改意見、審核論文；
郭猜：參與文獻調研及結果分析討論；
劉金凱：參與結果分析討論，協調辦理相關發表流程。