不同內(nèi)邊界條件下SEP事件實(shí)例的集合數(shù)值模擬及其在SEP 事件預(yù)報中的應(yīng)用

杜晨曦，敖先志，羅冰顯，王晶晶，李剛

（1.中國科學(xué)院國家空間科學(xué)中心，北京 100190；2.中國科學(xué)院大學(xué)，北京 100049；3.中國科學(xué)院空間態(tài)勢感知技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100190；4.美國阿拉巴馬大學(xué)漢茨維爾分校 空間科學(xué)系，漢茨維爾 AL35899）

0 引言

太陽高能粒子（solar energetic particles,SEPs）包括電子、質(zhì)子以及其他重離子，粒子能量范圍可以從幾十KeV 到GeV。太陽質(zhì)子事件（solar proton events,SPEs）屬于太陽高能粒子事件的一個子集，通常被定義為地球靜止軌道上＞10 MeV 的高能質(zhì)子通量在一段時間內(nèi)＞10 pfu。太陽質(zhì)子事件是太陽活動爆發(fā)所引起的災(zāi)害性空間天氣事件之一，其定量數(shù)值預(yù)報是空間環(huán)境態(tài)勢感知中非常重要的一個方面。

Forbush 報告了人類歷史上首次觀測到的太陽高能粒子事件[1]。太陽高能粒子事件通常被分為“脈沖型”（impulsive）和“漸進(jìn)型”（gradual）兩種類型[2]。兩種類型的SEP事件具有不同的觀測特征：“脈沖型”事件的通量曲線在觀測上體現(xiàn)出短時（幾個小時以內(nèi)）抬升的特征，而“漸進(jìn)型”事件的通量曲線則往往持續(xù)較長時間（2天甚至更長時間）。實(shí)際上，這兩種不同形態(tài)往往混在一起[3-5]，但總體來說大部分的SEP事件具有“漸進(jìn)型”特性[6]。一般認(rèn)為，“脈沖型”的SEP事件起源于耀斑爆發(fā)，而“漸進(jìn)型”的SEP事件起源于日冕物質(zhì)拋射（CME）所驅(qū)動的激波擴(kuò)散加速（diffusive shock acceleration,DSA）機(jī)制[7-9]。DSA 加速機(jī)制為一階費(fèi)米加速（first-order Fermiacceleration），由Axford、Krymskii、Blandford 和Bell 等在解釋銀河宇宙線的超新星激波加速過程時提出[10-15]。Drury對DSA 機(jī)制進(jìn)行了總結(jié)[16]。

Zank 等提出了一種基于DSA 加速機(jī)制的日球?qū)恿Ｗ蛹铀賯鬏斈Ｐ停碢ATH（the particle acceleration and transport in heliosphere）模型，并在隨后的數(shù)十年中對該模型進(jìn)行了不斷改進(jìn)和發(fā)展[17-25]。胡駿翔等對PATH 模型進(jìn)行了進(jìn)一步的改良和發(fā)展[26-27]，將原有模型中的激波演化從一維擴(kuò)展到二維，并建立了iPATH（improved PATH）模型。

本文在三個模型的基礎(chǔ)上，結(jié)合1 AU 處衛(wèi)星的太陽風(fēng)觀測參數(shù)和日冕儀的CME 觀測參數(shù)，建立了一套可用于預(yù)報SEP事件的數(shù)值方法。這三個模型包括：1）一個用來估算0.1 AU（大約20個太陽半徑(Rs)）處背景太陽風(fēng)參數(shù)的一維多方模型；2）反演CME 參數(shù)的冰激凌錐模型；3）iPATH 模型。利用該數(shù)值方法對發(fā)生于2013年9月30日的一次SEP實(shí)例進(jìn)行了數(shù)值模擬，并開展不同CME拋射速度和不同背景太陽風(fēng)溫度條件下的集合模擬試驗(yàn)。

1 模型簡介

1.1 iPATH模型簡述

iPATH 模型是對原有PATH 模型的進(jìn)一步發(fā)展，該數(shù)值模型使用了ZEUS-3D(v3.5)代碼[28-29]，在二維計算域中模擬太陽風(fēng)、CME 驅(qū)動的激波以及與之相關(guān)的粒子加速和傳輸過程。iPATH模型的計算分3個步驟：第1步，模擬背景太陽風(fēng)；第2步，在背景太陽風(fēng)中引入一個沖擊波來模擬CME并構(gòu)建了一個二維“洋蔥殼”（onion shell）以便計算粒子在激波前沿的加速過程；第3步，計算被激波加速的太陽高能粒子在行星際空間中的傳播過程。其中，第2步中所構(gòu)建的二維洋蔥殼及相應(yīng)的粒子加速計算過程是對原PATH 模型所進(jìn)行的主要改進(jìn)部分。

ZEUS-3D是一種基于FORTRAN語言編寫的非相對論磁流體力學(xué)（MHD）算子，采用交錯網(wǎng)格和迎風(fēng)格式，所求解的方程組為

式中：ρ是流體密度；V是流體速度；p是流體熱壓強(qiáng)；B是磁感應(yīng)強(qiáng)度；e是單位體積內(nèi)的能量，包括內(nèi)能和磁能；Φ是引力勢能。

ZEUS-3D本身具有在三維空間內(nèi)求解MHD方程組的能力，本文在iPATH 模型中第3個維度只取了有限個數(shù)的網(wǎng)格數(shù)量，即僅在一個平面內(nèi)（θ=90°：近似模擬赤道面或者黃道面）求解MHD方程組。這樣的處理方式雖然使得iPATH模型缺失了求解三維問題的能力，但是可以節(jié)省計算資源，減少計算時間，從而能夠在CME 爆發(fā)后在盡可能短的時間內(nèi)完成數(shù)值模擬計算，增加預(yù)警的時間提前量。

iPATH 模型以0.1 AU（約20Rs）為內(nèi)邊界，在內(nèi)邊界處引入一段時間的擾動（沖擊波）來模仿CME。內(nèi)邊界的背景太陽風(fēng)參數(shù)（太陽風(fēng)速度、磁場強(qiáng)度、密度和溫度）以及CME的參數(shù)（角寬度、拋射速度和擾動持續(xù)時間等）是驅(qū)動模型的必要參數(shù)。數(shù)值模擬的二維計算域，在0.1～2 AU 的徑向距離上均勻分布了1500個網(wǎng)格，在0°～360°的切向方向上均勻分布有360個網(wǎng)格。iPATH 模型計算時的初始磁場為如下帕克螺旋場，

式中：R0是源位置處的日心距；Br為磁場的徑向分量；Bφ為磁場的切向分量；B0為r=R0處的磁場；Ω為太陽自轉(zhuǎn)角速度；Usw為太陽風(fēng)速度；θ為球坐標(biāo)系中的仰角。

本文使用ACE 衛(wèi)星觀測的太陽風(fēng)參數(shù)，利用一維多方模型[30]估算內(nèi)邊界（0.1 AU）處的太陽風(fēng)參數(shù)并作為iPATH 模型的輸入來計算背景太陽風(fēng)。對于CME 擾動參數(shù)，我們使用了一套基于冰激凌錐模型的反演程序得到了內(nèi)邊界處CME拋射速度、角寬度、以及CME 源區(qū)位置的坐標(biāo)。

文獻(xiàn)[17-28]對粒子在激波附近的加速理論和數(shù)值模型已有非常詳細(xì)的討論。iPATH 模型相對于之前的PATH 模型而言，增加了對垂直擴(kuò)散系數(shù)的考慮，在激波波前上的不同位置都會有不同的擴(kuò)散系數(shù)，需要指出的是，此版本的iPATH 模型未考慮橫向擴(kuò)散對局地粒子分布函數(shù)的影響。被加速的高能粒子在激波附近發(fā)生逃逸，并沿著磁力線方向向外傳播，在傳播過程中受到行星際湍流影響而發(fā)生散射。粒子傳輸過程受到聚焦傳輸方程（focused transport equation）的控制，iPATH 模型的第三步使用了蒙特卡羅（Monte Carlo）方法來求解高能粒子的傳輸過程。

1.2 模型參數(shù)

發(fā)生于2013年9月30日的一次SEP事件被包括GOES15衛(wèi)星在內(nèi)的多顆空間科學(xué)衛(wèi)星所觀測到，事件表現(xiàn)為2天內(nèi)高能質(zhì)子通量的持續(xù)增高，并在第3天開始逐漸下降。此次事件是一次典型的由CME驅(qū)動激波加速粒子所引起的“漸進(jìn)型”SEP事件。9月29日22點(diǎn)12分（UT)），部署在太陽-地球引力平衡點(diǎn)L1處的SOHO衛(wèi)星觀測到了全暈CME事件的爆發(fā)。10月2日2時（UT）左右，ACE 衛(wèi)星觀測到CME 激波到達(dá)地球附近。GOES15衛(wèi)星所搭載EPEAD（Energetic Proton Electron and Alpha Detector）載荷的觀測顯示，在此次事件之前20多天的時間內(nèi)，高能粒子通量處于穩(wěn)定的低水平背景狀態(tài)，在事件結(jié)束后的20多天內(nèi)沒有觀測到典型的太陽高能粒子事件。這表明本文所研究的這一事件在時間上具有較好的孤立性，可以排除多個事件相互耦合所造成的復(fù)雜影響。SOHO衛(wèi)星的觀測顯示，在此次全暈CME事件發(fā)生的前后一段時間內(nèi)，沒有其他全暈CME 的發(fā)生。

如前所述，模型設(shè)置的網(wǎng)格為1500×360，為簡單起見，模型在計算背景太陽風(fēng)參數(shù)時假設(shè)內(nèi)邊界上所有點(diǎn)處的內(nèi)邊界條件都是軸對稱的，即同一物理量的數(shù)值在內(nèi)邊界所有的點(diǎn)上都是相同的。通過ACE衛(wèi)星的觀測獲得1 AU 附近（即地球附近）的背景太陽風(fēng)參數(shù)以及采用一維多方模型計算得到的內(nèi)邊界處的背景太陽風(fēng)參數(shù)分別如表1所示。其中在計算磁場分量時，假設(shè)磁場在z方向的分量很小，因此磁場分量的計算結(jié)果簡化為二維平面中的2個分量，即徑向分量和切向分量。

表1 背景太陽風(fēng)參數(shù)Table 1 Background solar wind parameters

本文利用冰激凌錐模型[31]來反演CME參數(shù)。首先使用人工方法對SOHO衛(wèi)星日冕儀的白光日冕圖像進(jìn)行CME 前沿的信息提取。圖1分別展示了對不同時刻的C3圖像進(jìn)行人工CME信息提取的過程，其中實(shí)心點(diǎn)為人工標(biāo)記的CME 前沿。

圖1 人工CME 前沿標(biāo)記Fig.1 Artificial CME frontier marking

進(jìn)而通過計算，由反演程序輸出了圖2所示的擬合曲線，并獲得了表2中的CME 參數(shù)。圖2顯示了冰激凌錐模型的擬合曲線、CME 拋射速度和CME 拋射的角寬度。

圖2 冰激凌錐模型擬合曲線Fig.2 Fitting result of the ice-cream cone model

表2 冰激凌錐模型反演結(jié)果Table 2 Inversion result of the ice-cream cone model

iPATH 模型中的粒子注入率為[25]

2 數(shù)值模擬結(jié)果

2.1 iPATH模型實(shí)例模擬結(jié)果

利用上述方法得到的模型輸入?yún)?shù)，以地球位置為觀測點(diǎn)進(jìn)行了數(shù)值模擬計算。在模型中設(shè)定CME 中心方位角為90°，并通過冰激凌錐模型反演的源區(qū)位置計算了地球在坐標(biāo)系中相應(yīng)的位置。圖3為歸一化數(shù)密度的數(shù)值模擬結(jié)果，分別展示了4 個不同時刻的CME 和激波在黃道面內(nèi)傳播的過程中磁力線彎曲的狀態(tài)和演化。圖3 中：綠色標(biāo)記點(diǎn)為地球位置；螺旋實(shí)線為磁力線，其中綠色實(shí)線為穿過地球的磁力線；顏色棒代表了歸一化的密度值，顏色越深代表數(shù)值越大。由圖可見，觀測點(diǎn)與激波波前處磁力線連接的位置對太陽高能粒子的模擬結(jié)果有著非常重要的影響；激波在傳播的過程中，逐漸由準(zhǔn)平行激波轉(zhuǎn)化為準(zhǔn)垂直激波，同時連接點(diǎn)的位置也從激波邊緣逐漸向激波中心移動。從圖3(c)可以看出，激波在10月2日13：42：45時刻已經(jīng)越過了地球所在位置，這與ACE 所觀測到的激波到達(dá)時間基本相符。

圖3 CME 和激波在黃道面中的傳播Fig.3 The simulated propagation of the CME and the shock in the ecliptic plane

圖4對比了此次太陽高能粒子事件的數(shù)值模擬結(jié)果與GOES15衛(wèi)星的觀測結(jié)果，其中：紅色、藍(lán)色、黑色的曲線分別代表能量≥10 MeV、≥50 MeV和≥100 MeV 高能質(zhì)子的積分通量；實(shí)線和虛線分別代表GOES15的觀測結(jié)果和數(shù)值模擬的結(jié)果。由于CME 爆發(fā)在太陽表面的附近，而模型的內(nèi)邊界在0.1 AU 處，所以模型的計算結(jié)果忽略了CME從太陽表面附近傳播到0.1 AU 處的過程。在對比時，通過反演得到的CME 拋射速度估算了CME傳播到0.1 AU 處的時間，從而對模擬結(jié)果中高能粒子事件的開始時間進(jìn)行了簡單的校正。由圖可見：數(shù)值模擬結(jié)果中＞10 MeV 高能粒子通量與觀測結(jié)果吻合較好，其變化趨勢與觀測結(jié)果較為相近；對于能量＞50 MeV 和＞100 MeV 高能質(zhì)子的積分通量而言，模擬結(jié)果雖然能夠大致體現(xiàn)出與觀測結(jié)果相應(yīng)的變化趨勢，但其通量高于觀測值，尤其是能量＞100 MeV 的高能質(zhì)子積分通量差異顯著。造成這種結(jié)果的可能原因之一是，由于在模擬中對粒子注入率的處理與實(shí)際物理過程有很大差異，導(dǎo)致較高能段的粒子注入率偏高。另外一種可能的因素是二維數(shù)值模擬不能很好地模擬三維空間中的CME傳播，導(dǎo)致激波壓縮比計算值偏高。

圖4 GOES衛(wèi)星觀測的SEP 事件單向積分通量與模擬結(jié)果對比Fig.4 Comparison between the observed unidirectional integrated flux obtained by GOES satellite and the simulated results

2.2 不同CME 拋射速度條件下的集合模擬結(jié)果

限于觀測手段，錐模型擬合的CME 的拋射速度和主拋射方向往往和實(shí)際情況相差較大。為了研究不同的CME拋射速度對iPATH模型數(shù)值模擬結(jié)果的影響，本文分別選取了CME 拋射速度分別為900、1000、1100和1300 km/s進(jìn)行數(shù)值模擬，并將模擬結(jié)果與前述模擬結(jié)果進(jìn)行對比。圖5顯示了同一時刻不同CME拋射速度條件下的CME/激波狀態(tài)，從圖5(a)到圖5(d)，CME拋射速度依次增大，其中圖5(c)為2.1節(jié)中的模擬結(jié)果，與圖3(c)相同?？梢钥吹?，在同一時刻CME拋射速度越高，激波波前距離太陽越遠(yuǎn)。圖6顯示了不同CME拋射速度條件下能量＞10 MeV 高能質(zhì)子的單向積分通量?？梢钥吹?，積分通量曲線在不同CME拋射速度下有較為明顯的差異：在事件開始的前段，積分通量呈現(xiàn)出隨著CME 拋射速度的增大而增高的趨勢；在事件后段，積分通量的差異開始收縮。

圖5 不同CME 拋射速度條件下的CME/激波對比Fig.5 Comparison of CME/shock propagations of different initial CME speeds

圖6 不同CME 拋射速度條件下＞10 MeV 的高能質(zhì)子通量模擬結(jié)果對比Fig.6 Comparison of simulated results of energetic flux of particles of more than 10 MeV with different initial CME speeds

進(jìn)一步考查不同CME 拋射速度條件下高能粒子能譜的特征。圖7中所示為不同CME拋射速度條件下整個事件的高能粒子能譜分布以及與之相對應(yīng)的分段線性擬合的結(jié)果，其中：不同顏色的實(shí)心點(diǎn)代表不同拋射速度，直線代表對能譜進(jìn)行分段擬合的結(jié)果?？梢钥闯觯珻ME拋射速度越大的高能粒子能譜越高。

圖7 不同CME 拋射速度條件下高能質(zhì)子全事件積分能譜分布及擬合結(jié)果Fig.7 Event-integrated power spectra and fitting results for different initial ejection speeds

為了定量顯示高能粒子能譜間的差異，使用了形如y=kx+b的函數(shù)對高能粒子能譜進(jìn)行了最小二乘擬合，并給出了與之相應(yīng)的分段擬合結(jié)果，如表3所示。表中給出了擬合直線的k和b，其中下角標(biāo)1、2分別表示第1段和第2段的擬合直線。結(jié)果顯示粒子能譜間具有較大差異。

表3 不同CME拋射速度條件粒子能譜擬合結(jié)果Table 3 Fitting result of the particle spectrum for different CME ejection speeds

2.3 不同內(nèi)邊界背景太陽風(fēng)溫度條件下的集合模擬結(jié)果

一般情況下，內(nèi)日冕的太陽風(fēng)溫度無法直接測量。為進(jìn)一步研究不同內(nèi)邊界背景太陽風(fēng)溫度對iPATH 模型數(shù)值模擬結(jié)果的影響，采取與2.2節(jié)相似的研究方法。在此次事件中，我們通過一維多方模型估算的內(nèi)邊界背景太陽風(fēng)溫度為0.4287 MK。本文分別選取此數(shù)值的0.5倍、0.75倍、1.5倍和2 倍作為內(nèi)邊界溫度條件進(jìn)行數(shù)值模擬，并將結(jié)果與前述模擬結(jié)果進(jìn)行對比。圖8所示為同一時刻不同內(nèi)邊界背景太陽風(fēng)溫度條件下的CME 激波狀態(tài)，由圖8(a)到圖8(d)溫度依次增大，其中圖8(c)為前述實(shí)例模擬的結(jié)果?？梢钥吹?，在同一時刻背景太陽風(fēng)溫度越高，激波波前面距離太陽越遠(yuǎn)，但這種變化相對于圖5來說明顯更小。這說明背景太陽風(fēng)溫度的增高會引起CME 速度的微弱變化。圖9 所示為不同內(nèi)邊界背景太陽風(fēng)溫度條件下能量＞10 MeV高能質(zhì)子的單向積分通量。可以看到：積分通量曲線在不同內(nèi)邊界背景太陽風(fēng)溫度下的差異很??；在事件開始的前段，積分通量呈現(xiàn)出隨著內(nèi)邊界背景太陽風(fēng)溫度增大而有微弱增高的趨勢，在事件后段則呈現(xiàn)出與前段相反的趨勢。

進(jìn)一步考查不同內(nèi)邊界背景太陽風(fēng)溫度條件下高能粒子能譜的特征。圖10所示為不同內(nèi)邊界背景太陽風(fēng)溫度條件下的高能粒子能譜分布以及與之相對應(yīng)的分段線性擬合的結(jié)果，其中：不同顏色的實(shí)心點(diǎn)代表內(nèi)邊界處不同的背景太陽風(fēng)溫度，直線代表對能譜進(jìn)行分段擬合的結(jié)果。可以看出，內(nèi)邊界背景太陽風(fēng)溫度的變化對高能粒子能譜的影響很小。為了定量地顯示高能粒子能譜間的差異，本文仍然使用形如y=kx+b的函數(shù)對高能粒子能譜進(jìn)行了最小二乘擬合，并給出了與之相應(yīng)的分段擬合結(jié)果，如表4所示。結(jié)果清晰地表明內(nèi)邊界背景太陽風(fēng)溫度的變化對SEP事件中高能粒子能譜的影響很小。

圖8 內(nèi)邊界不同背景太陽風(fēng)溫度條件下的CME/激波對比Fig.8 Comparison of CME/shock propagation at different background solar wind temperatures

圖9 內(nèi)邊界不同背景太陽風(fēng)溫度條件下＞10 MeV 的高能質(zhì)子通量模擬結(jié)果對比Fig.9 Comparison of the simulated results of energetic flux of particle of more than 10 MeV at different background solar wind temperatures

圖10 內(nèi)邊界不同背景太陽風(fēng)溫度條件下的高能質(zhì)子全事件積分能譜分布及擬合結(jié)果Fig.10 Event-integrated power spectra and fitting results for different background solar wind temperatures

表4 不同背景太陽風(fēng)溫度條件下粒子能譜擬合結(jié)果Table 4 Fitting result of the particle spectrum at different background solar wind temperatures

3 結(jié)束語

本文使用iPATH 模型對2013年9月30日的一次孤立SEP事件進(jìn)行了數(shù)值模擬。數(shù)值模擬結(jié)果和GOES衛(wèi)星觀測結(jié)果的對比顯示：數(shù)值模擬的能量＞10 MeV 的高能粒子的通量和觀測值吻合較好，而能量＞100 MeV 的高能粒子通量的模擬結(jié)果高于觀測值。造成這種情況的可能原因之一是不同能量的高能粒子在加速過程中經(jīng)歷不一樣的粒子注入過程，或者說注入率不一樣。另外一種可能的原因是我們在計算過程中采用的是二維模型，因此對激波參數(shù)的計算和實(shí)際的三維情況存在偏差。另一方面，基于SOHO衛(wèi)星白光日冕觀測圖像，通過冰激凌錐模型反演CME參數(shù)，然而模型反演的CME 參數(shù)往往和實(shí)際情況相差較大。鑒于SEP事件的形成過程受到CME 初始拋射速度的影響較大，為了定量地研究CME 初始拋射速度帶來的影響，進(jìn)一步開展了不同CME 初始拋射速度下的集合模擬試驗(yàn)。集合模擬的結(jié)果表明更高的CME 初始拋射速度會引起更強(qiáng)的SEP事件。

我們的預(yù)報模型采用了一個較為簡單的一維多方模型來估計內(nèi)邊界處背景太陽風(fēng)的參數(shù)。一般而言，內(nèi)邊界處背景太陽風(fēng)的速度、密度和磁場的估算較為容易，而溫度的估算則相對困難，一維多方模型所估算的溫度往往高于實(shí)際的溫度。此外，以往的文獻(xiàn)基本上沒有關(guān)于背景太陽風(fēng)溫度對于SEP事件影響的相關(guān)討論。因此，我們開展了不同背景太陽風(fēng)溫度條件下的集合模擬研究。數(shù)值計算的結(jié)果表明：內(nèi)邊界背景太陽風(fēng)溫度的變化對于地球附近SEP事件中高能粒子通量和能譜的影響幾乎可以忽略不計。該結(jié)論在空間環(huán)境態(tài)勢感知中具有重要的應(yīng)用意義。在為各類空間探測計劃提供保障服務(wù)、實(shí)際進(jìn)行太陽質(zhì)子事件預(yù)報的過程中，即使我們不能非常準(zhǔn)確地獲取靠近太陽的內(nèi)邊界處背景太陽風(fēng)的溫度，也不會顯著地影響預(yù)報的結(jié)果。

致謝

本文使用了ICA(www.ica.smu.ca)的CLARKE在NSERC的支持下開發(fā)的ZEUS-3D模型，在此表示感謝。（Use of ZEUS-3D,developed by CLARKE D.A.at the ICA(www.ica.smu.ca)with support from NSERC,isacknowledged.）