氣象擴散耦合模型對經典示蹤試驗ETEX-1的模擬分析

郭 歡,張俊芳,李鳴野,黃 莎,趙多新

(中國輻射防護研究院 核環(huán)境模擬與評價技術重點實驗室,太原 030006)

核事故發(fā)生后,氣載放射性污染物在大氣中迅速擴散,短時間內可能造成大范圍的環(huán)境影響,若能精確預測污染物的擴散規(guī)律,就可為及時制定應急響應決策提供參考[1]。在評估污染物在大氣中的擴散過程時,單純依靠觀測資料存在數(shù)據(jù)有限、無法預報未來擴散趨勢等短板問題,而數(shù)值模擬因其不受試驗條件和天氣條件限制等優(yōu)點成為開展大氣擴散研究的一種常用手段[2]。大氣擴散主要受大氣邊界層中的湍流控制,而湍流存在隨機性,通常采用統(tǒng)計學方法將野外示蹤試驗數(shù)據(jù)與數(shù)值模式的模擬結果作比對來描述其對擴散的影響[3-4]。陳曉秋[5]利用大亞灣核電廠的示蹤試驗數(shù)據(jù)對隨機游走模式的數(shù)值模擬能力進行了驗證;胡二邦等人[6]基于福建惠安核電廠址的野外示蹤試驗數(shù)據(jù)對高斯煙羽模式在復雜濱海地形下的模擬進行了分析研究;王博[7]以野外示蹤試驗數(shù)據(jù)為參考,對CFD的大氣擴散特征模擬能力進行了驗證;朱好等人[8]利用湖南某廠址的野外示蹤試驗數(shù)據(jù)對拉格朗日煙團模型CALPUFF在復雜地形下的近場適用性開展了研究。

本文采用WRF V3.5中尺度氣象模式和FLEXPART-WRF V3.1粒子擴散模式,并根據(jù)WRF版本,調節(jié)WRF輸出數(shù)據(jù)中的默認氣象要素以及調節(jié)模式輸出時空分辨率、模擬區(qū)域、物理參數(shù)化方案等參數(shù)信息,改動接口程序以適應V3.1版本粒子擴散模式的輸入需求,利用國外經典示蹤試驗ETEX-1數(shù)據(jù)集和由丹麥國家環(huán)境研究所(NERI)開發(fā)的MVK工具中的統(tǒng)計方法,對試驗期間的情景進行了數(shù)值再現(xiàn),并將模擬結果與相關的氣象觀測數(shù)據(jù)及現(xiàn)場濃度數(shù)據(jù)作比對,對耦合模式模擬效果進行了驗證分析。

1 模式設置與研究方法

1.1 ETEX試驗概況

歐洲長距離示蹤試驗(ETEX)[9]包括在1994年10月和11月期間開展的兩次試驗,根據(jù)氣象背景統(tǒng)計分析,中歐上空出現(xiàn)西風氣流時大多數(shù)歐洲國家會受到示蹤羽流的影響,天氣條件最有利于試驗開展,因此釋放點定于法國西部的蒙特菲爾。為了避免交叉污染和樣本污染,第一次釋放PMCH(全氟甲基環(huán)己烷)作為示蹤物,第二次釋放PMCP(全氟甲基環(huán)戊烷),且兩次試驗的時間間隔超過2周。試驗中煙羽遷移范圍達2 000余km,在釋放后72小時內每間隔3小時進行采樣,使用遍布歐洲17個國家的168個地面監(jiān)測站點的采樣網(wǎng)收集空氣樣本,兩次實驗采集到約9 000份樣本。示蹤數(shù)據(jù)庫中匯集了釋放特征信息、地面和高空的常規(guī)與附加氣象觀測、地面和高空濃度等數(shù)據(jù)。第一次釋放(ETEX-1)期間,天氣狀態(tài)平穩(wěn),沒有出現(xiàn)高壓或低壓中心,也沒有延伸的脊或槽靠近釋放點,無鋒面系統(tǒng)經過釋放地點,在穩(wěn)定的強西風和西南風氣流引導下,示蹤羽流被輸送到整個歐洲。第二次試驗(ETEX-2)期間受釋放初期鋒面系統(tǒng)過境的影響,天氣條件惡劣,煙羽并未按照預期覆蓋所布監(jiān)測點,對試驗的開展和空氣樣本數(shù)據(jù)質量都產生了極大不利影響,國際上通常認為第一次試驗較為成功,因此本文選擇第一次試驗(ETEX-1)期間的情形作模擬分析。

圖1為示蹤試驗期間的站點分布,包括歐洲168站和法國西北部81站,其中,歐洲168站既進行氣象觀測又收集濃度樣本,法國81站僅進行常規(guī)氣象觀測。

紅:釋放點;藍:歐洲168站;綠:法國81站。

1.2 WRF 模式參數(shù)設置

根據(jù)示蹤試驗時段,選用歐洲中期天氣預報中心(European Centre for Medium-Range Weather Forecasts,ECMWF)的再分析數(shù)據(jù)作為WRF模式的初始場和邊界條件,對試驗期間的72 h進行模擬,即1994年10月23日12時至1994年10月26日12時(UTC),數(shù)據(jù)空間分辨率為1°×1°,地形數(shù)據(jù)采用美國地理學會(USGS)的10′分辨率地形數(shù)據(jù)。對模擬的歐洲地區(qū)進行網(wǎng)格劃分,模擬中心位置設為49°N,12°E,采用兩層嵌套,第一層(d01)網(wǎng)格分辨設定為30 km,模擬區(qū)網(wǎng)格數(shù)為111×111,第二層(d02)網(wǎng)格分辨設定為10 km,模擬區(qū)網(wǎng)格數(shù)為231×231。在對物理參數(shù)化方案的設置上,微物理參數(shù)化方案選擇WSM3,長波和短波輻射方案分別選擇RRTM和Dudhia,邊界層參數(shù)化方案選擇MYNN2.5,陸面過程方案選NOAH。

1.3 FLEXPART-WRF模式參數(shù)設置

FLEXPART-WRF由驅動場氣象模型WRF和主體擴散模型FLEXPART組成。WRF模型以原始全球數(shù)值預報產品作為輸入,經過模式動力降尺度后可獲取更高時空分辨率的高質量氣象場,從而為擴散模型提供更精確的氣象初始場。FLEXPART模式由挪威大氣研究所開發(fā),是一種典型的拉格朗日粒子擴散模式,通過設置釋放大量隨機粒子,跟蹤其在湍流中的運動來獲取粒子位置和擴散濃度[10]。FLEXPART在模擬過程中引入邊界層參數(shù)化方案來反映對流和湍流運動特征,因此可以更準確地描述粒子在大氣中的運動過程。兩種模式經過耦合,通過改善氣象場模擬效果可提高擴散模擬精度。但耦合模型也存在一定的局限性,其描述的是幾百至幾千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時間),連續(xù)釋放12小時。

1.4 驗證方法

NERI以Hanna等[11]提出的定量統(tǒng)計參數(shù)為基礎,研發(fā)了一套模型驗證比對工具MVK(Model Validation Kit),其中的定量分析程序(Boot)包含全面的ASTM(American Society for Testing and Materials)模型驗證評價方法,通過將試驗觀測值和模型預測值進行定量統(tǒng)計分析,可以有效評價模型的準確性。統(tǒng)計指標包括平均值(Mean)、平均偏差(Bias)、標準偏差(Sigma)、偏差分數(shù)(FB)、幾何平均偏差(MG)、幾何平均方差(VG)、標準化均方誤差(NMSE)、相關系數(shù)(R)等。另外,美國國家海洋和大氣管理局(NOAA)組織的DATEM項目專門用于大尺度大氣擴散模式的有效性評估,評估中使用相關系數(shù)(R)、偏差分數(shù)(FB)、空間吻合度系數(shù)(FMS)、KS檢驗(Kolmogorov-Smirnov,KS)組合為一個模式評分指標RANK對模式有效性進行評價。本文結合MVK工具和NOAA對模式的評分指標評價方法,對耦合模式進行了ETEX-1試驗的再現(xiàn)模擬分析,并將分析結果與美國空氣資源實驗室(ARL)官網(wǎng)公布的HYSPLIT模式對ETEX-1模擬結果的RANK指標作對比。相關計算公式如表2所示。

表2 指標說明

2 結果分析

2.1 風場模擬結果分析

在核事故應急預警中,濃度場是重要的評估對象,而風場作為大氣擴散的驅動因子會直接影響濃度分布,因此首先對風場進行分析。圖2表示不同時刻模擬的地面風場。ETEX-1試驗覆蓋面廣,范圍涉及歐洲多個國家和地區(qū)。與示蹤試驗相適應,模擬區(qū)域覆蓋整個西歐以及部分東歐地區(qū),且包含北大西洋的一部分、北海、地中海、波羅的海等海域,試驗期間模擬區(qū)西部有來自北大西洋的海風,在地中海、北海附近多出現(xiàn)環(huán)流。在試驗過程中的大多數(shù)時刻,釋放點均處于北大西洋吹向陸面的海風之下,且受地形和背景天氣形勢的共同影響,氣流在法國東部開始分流,一部分融入北海環(huán)流中,一部分流向地中海。定性分析結果認為風場較好地模擬了地面風場受海陸影響而產生的擾動變化情況。

圖2 不同時刻的模擬地面風場

對歐洲168個常規(guī)地面氣象站、法國西北部81個常規(guī)地面氣象站的觀測結果與模擬的地面風場分別進行對比,對比結果見圖3～圖4,法國站U、V風矢量的模擬值與觀測值之間的相關性比歐洲站高,均方根誤差比歐洲站低,除個別時刻外,相關系數(shù)和均方根誤差隨時間的波動幅度較小。統(tǒng)計分析全部站點數(shù)據(jù)后,得出U、V風分量的相關系數(shù)均值為0.43,均方根誤差均值為3.4 m/s,總體上,風場模擬結果能夠代表實際風場情況。由于擴散與風場直接相關,因此風場從一定程度上為擴散模擬奠定了良好的氣象條件基礎。

圖3 U、V風分量相關系數(shù)隨時間的變化

圖4 U、V風分量均方根誤差隨時間的變化

2.2 擴散場模擬結果分析

擴散與風場形勢密切相關,示蹤劑被釋放后,隨著模擬時間的推移,示蹤粒子跟隨風場向遠距離遷移擴散,且示蹤物分布范圍逐漸變廣。圖5(a)給出示蹤粒子擴散所造成的近地面濃度分布情況,結合模擬的地面風場分布進行分析,從釋放點到遠距離的擴散過程中,粒子呈現(xiàn)主要往東北和部分向東南遷移的趨勢,近地面過境濃度分布與模擬風場結果匹配較好,從煙羽覆蓋區(qū)和采樣點的分布上也可以看出,模式對煙羽的空間位置捕捉效果較好。圖5(b)為示蹤粒子在整個試驗階段的最大過境濃度分布,其向東北的走勢與風場形勢吻合度很高,且煙云覆蓋區(qū)域與歐洲相應的最大地面濃度監(jiān)測站點的分布近乎一致。此外,根據(jù)模式計算結果,模擬最大濃度出現(xiàn)在釋放點附近,與實際相符,說明耦合模式模擬結果能較好地反映風場對氣載污染物的聚集和擴散的影響。

圖5 近地面濃度分布

統(tǒng)計濃度模擬值與168站濃度監(jiān)測值之間的平均絕對偏差,以及各自濃度不為0的站點個數(shù),如表3所示。其中,時長表示從釋放開始每隔3 h排序,這是為了與3 h時間分辨的試驗采樣數(shù)據(jù)保持一致。結果表明:在前24 h,觀測與模擬的平均絕對偏差呈現(xiàn)先增大后減小的趨勢,最高偏差為1.37 ng/m3,出現(xiàn)在釋放后的第15 h。第27 h至第36 h出現(xiàn)短暫偏差增大的過程,隨后逐漸較小,趨于平穩(wěn)。對于濃度不為0的統(tǒng)計結果,觀測和模擬個數(shù)均逐漸增大,這是由于示蹤粒子一直在運動遷移,必然逐漸覆蓋新的區(qū)域。結合模擬濃度的空間分布情況,釋放63 h后,煙羽主要部分超出了采樣網(wǎng)捕捉范圍,因此,中歐雖仍存在低濃度測量值,但模擬濃度不為0的個數(shù)相比觀測明顯減少。

表3 濃度平均絕對偏差和濃度不為0的點位數(shù)據(jù)統(tǒng)計

散點圖和分位數(shù)圖在一定程度上能反映模型預測情況。在散點圖中,觀測值和預測值成對出現(xiàn),可以直觀顯示模型預測過高或過低的程度。在統(tǒng)計分析時,剔除觀測數(shù)據(jù)中標記為“-0.88”和“-0.99”的數(shù)據(jù),其中,“-0.88”表示該站點位置能夠檢測到濃度,但未能給出有效數(shù)據(jù),“-0.99”表示該站未檢測到濃度值。圖6中,數(shù)據(jù)點多落于左上角區(qū)域,說明模式模擬值比觀測值整體偏高。分位數(shù)圖(圖7)將每個觀測和預測數(shù)據(jù)從最低到最高進行排序,可以反映模式測預測趨勢,從圖中也可看出,模式對濃度的模擬總體偏高。

圖6 模型散點圖

圖7 模型分位數(shù)圖

選取2個指標分析其隨時間的變化,從時間分布的角度考察模式性能(如圖8所示)。Bias表征預測均值和觀測均值的差值,從圖8來看,除個別時刻Bias為負值,在模擬的絕大多數(shù)時刻均為正值,到后期二者差距逐漸接近于0,說明模式對濃度的模擬整體偏高。NMSE在釋放開始后第3 h至第18 h較大,這可能是由于在短時間內釋放粒子數(shù)較多,模擬擴散不充分造成局部濃度積聚,從而導致觀測與模擬間產生了較大差距,之后減小逐漸穩(wěn)定趨于0。

圖8 模擬濃度與觀測值的Bias、NMSE隨時間的分布

計算統(tǒng)計學指標匯總如表4所示。FB和MG均可反映系統(tǒng)誤差,FB為無量綱數(shù),常用于對多量級濃度的比較,個別觀測或預測濃度出現(xiàn)極大值會對其產生較大影響,因而在觀測和預測完全相異時FB仍有可能為0,而MG因采用對數(shù)法對高值不敏感,但觀測或模擬值近0會對其產生顯著影響。VG反映由系統(tǒng)誤差導致的幾何方差和非系統(tǒng)誤差導致的幾何方差之和,對于預測濃度和觀測濃度多量級變化的情況,更宜評估MG和VG。NMSE為偏差統(tǒng)計參量,對觀測和預測間的偏差表現(xiàn)敏感。理想狀態(tài)下,數(shù)值模擬的MG、VG值為1,FB和NMSE為0,但由于自然環(huán)境中大氣運動存在隨機性,預測結果很難達到此標準。FMS表示濃度觀測區(qū)域和預測區(qū)域重疊的百分比,反映了模式的空間模擬效果[12]。KS(Kolomogorov-Smirnov)定義為兩個累積分布之間的最大差值,是唯一根據(jù)未配對數(shù)據(jù)計算的統(tǒng)計參量,可衡量模型對實際觀測濃度分布的再現(xiàn)程度,任何兩個分布之間的最大差異不能超過100%。Chang等[13]綜合考慮試驗情景復雜程度以及氣象條件的隨機波動等影響因素,在不要求模擬值和實測值在時間和空間一一對應的條件下,推薦了指標的可接受范圍,即-0.3

3 結論與討論

本文將中尺度氣象模型WRF和拉格朗日粒子擴散模型FLEXPART-WRF適應性耦合后進行了ETEX-1試驗的再現(xiàn)數(shù)值模擬,以試驗數(shù)據(jù)集為參考,結合MVK統(tǒng)計學指標,將模擬結果與氣象觀測和濃度監(jiān)測數(shù)據(jù)分別作對比分析。從擴散整體情況看,粒子擴散趨勢與流場相匹配,近地面過境濃度覆蓋大部分監(jiān)測站點;模擬與觀測之間的偏差、標準化均方誤差等指標值基本符合國際標準,說明模式有效地模擬了試驗期間的天氣過程和區(qū)域流場及擴散特征。本次模擬得到RANK值為1.89,比ARL實驗室官網(wǎng)HYSPLIT對ETEX-1模擬結果的RANK小0.07,二者非常接近,認為氣象擴散模型相耦合可有效反映放射性物質的遷移擴散規(guī)律,可用于輔助評價核事故下氣載放射性污染物擴散對環(huán)境的影響。

需要說明的是,為了適應對ETEX-1試驗的再現(xiàn)需求,本文采用再分析氣象數(shù)據(jù)作為初始數(shù)據(jù),在實際應急時,模式也可利用實時預報數(shù)據(jù)開展預報分析。此外,在未來的研究工作中,還需重點關注WRF模式中物理參數(shù)化過程的優(yōu)化分析,以及采用多種氣象觀測資料進行數(shù)據(jù)同化以期獲取更精確的氣象場,從而進一步改善擴散模式的模擬效果。