不同内边界条件下SEP事件实例的集合数值模拟及其在SEP 事件预报中的应用

杜晨曦，敖先志，罗冰显，王晶晶，李刚

（1.中国科学院国家空间科学中心，北京 100190；2.中国科学院大学，北京 100049；3.中国科学院空间态势感知技术重点实验室，北京 100190；4.美国阿拉巴马大学汉茨维尔分校 空间科学系，汉茨维尔 AL35899）

0 引言

太阳高能粒子（solar energetic particles,SEPs）包括电子、质子以及其他重离子，粒子能量范围可以从几十KeV 到GeV。太阳质子事件（solar proton events,SPEs）属于太阳高能粒子事件的一个子集，通常被定义为地球静止轨道上＞10 MeV 的高能质子通量在一段时间内＞10 pfu。太阳质子事件是太阳活动爆发所引起的灾害性空间天气事件之一，其定量数值预报是空间环境态势感知中非常重要的一个方面。

Forbush 报告了人类历史上首次观测到的太阳高能粒子事件[1]。太阳高能粒子事件通常被分为“脉冲型”（impulsive）和“渐进型”（gradual）两种类型[2]。两种类型的SEP事件具有不同的观测特征：“脉冲型”事件的通量曲线在观测上体现出短时（几个小时以内）抬升的特征，而“渐进型”事件的通量曲线则往往持续较长时间（2天甚至更长时间）。实际上，这两种不同形态往往混在一起[3-5]，但总体来说大部分的SEP事件具有“渐进型”特性[6]。一般认为，“脉冲型”的SEP事件起源于耀斑爆发，而“渐进型”的SEP事件起源于日冕物质抛射（CME）所驱动的激波扩散加速（diffusive shock acceleration,DSA）机制[7-9]。DSA 加速机制为一阶费米加速（first-order Fermiacceleration），由Axford、Krymskii、Blandford 和Bell 等在解释银河宇宙线的超新星激波加速过程时提出[10-15]。Drury对DSA 机制进行了总结[16]。

Zank 等提出了一种基于DSA 加速机制的日球层粒子加速传输模型，即PATH（the particle acceleration and transport in heliosphere）模型，并在随后的数十年中对该模型进行了不断改进和发展[17-25]。胡骏翔等对PATH 模型进行了进一步的改良和发展[26-27]，将原有模型中的激波演化从一维扩展到二维，并建立了iPATH（improved PATH）模型。

本文在三个模型的基础上，结合1 AU 处卫星的太阳风观测参数和日冕仪的CME 观测参数，建立了一套可用于预报SEP事件的数值方法。这三个模型包括：1）一个用来估算0.1 AU（大约20个太阳半径(Rs)）处背景太阳风参数的一维多方模型；2）反演CME 参数的冰激凌锥模型；3）iPATH 模型。利用该数值方法对发生于2013年9月30日的一次SEP实例进行了数值模拟，并开展不同CME抛射速度和不同背景太阳风温度条件下的集合模拟试验。

1 模型简介

1.1 iPATH模型简述

iPATH 模型是对原有PATH 模型的进一步发展，该数值模型使用了ZEUS-3D(v3.5)代码[28-29]，在二维计算域中模拟太阳风、CME 驱动的激波以及与之相关的粒子加速和传输过程。iPATH模型的计算分3个步骤：第1步，模拟背景太阳风；第2步，在背景太阳风中引入一个冲击波来模拟CME并构建了一个二维“洋葱壳”（onion shell）以便计算粒子在激波前沿的加速过程；第3步，计算被激波加速的太阳高能粒子在行星际空间中的传播过程。其中，第2步中所构建的二维洋葱壳及相应的粒子加速计算过程是对原PATH 模型所进行的主要改进部分。

ZEUS-3D是一种基于FORTRAN语言编写的非相对论磁流体力学（MHD）算子，采用交错网格和迎风格式，所求解的方程组为

式中：ρ是流体密度；V是流体速度；p是流体热压强；B是磁感应强度；e是单位体积内的能量，包括内能和磁能；Φ是引力势能。

ZEUS-3D本身具有在三维空间内求解MHD方程组的能力，本文在iPATH 模型中第3个维度只取了有限个数的网格数量，即仅在一个平面内（θ=90°：近似模拟赤道面或者黄道面）求解MHD方程组。这样的处理方式虽然使得iPATH模型缺失了求解三维问题的能力，但是可以节省计算资源，减少计算时间，从而能够在CME 爆发后在尽可能短的时间内完成数值模拟计算，增加预警的时间提前量。

iPATH 模型以0.1 AU（约20Rs）为内边界，在内边界处引入一段时间的扰动（冲击波）来模仿CME。内边界的背景太阳风参数（太阳风速度、磁场强度、密度和温度）以及CME的参数（角宽度、抛射速度和扰动持续时间等）是驱动模型的必要参数。数值模拟的二维计算域，在0.1～2 AU 的径向距离上均匀分布了1500个网格，在0°～360°的切向方向上均匀分布有360个网格。iPATH 模型计算时的初始磁场为如下帕克螺旋场，

式中：R0是源位置处的日心距；Br为磁场的径向分量；Bφ为磁场的切向分量；B0为r=R0处的磁场；Ω为太阳自转角速度；Usw为太阳风速度；θ为球坐标系中的仰角。

本文使用ACE 卫星观测的太阳风参数，利用一维多方模型[30]估算内边界（0.1 AU）处的太阳风参数并作为iPATH 模型的输入来计算背景太阳风。对于CME 扰动参数，我们使用了一套基于冰激凌锥模型的反演程序得到了内边界处CME抛射速度、角宽度、以及CME 源区位置的坐标。

文献[17-28]对粒子在激波附近的加速理论和数值模型已有非常详细的讨论。iPATH 模型相对于之前的PATH 模型而言，增加了对垂直扩散系数的考虑，在激波波前上的不同位置都会有不同的扩散系数，需要指出的是，此版本的iPATH 模型未考虑横向扩散对局地粒子分布函数的影响。被加速的高能粒子在激波附近发生逃逸，并沿着磁力线方向向外传播，在传播过程中受到行星际湍流影响而发生散射。粒子传输过程受到聚焦传输方程（focused transport equation）的控制，iPATH 模型的第三步使用了蒙特卡罗（Monte Carlo）方法来求解高能粒子的传输过程。

1.2 模型参数

发生于2013年9月30日的一次SEP事件被包括GOES15卫星在内的多颗空间科学卫星所观测到，事件表现为2天内高能质子通量的持续增高，并在第3天开始逐渐下降。此次事件是一次典型的由CME驱动激波加速粒子所引起的“渐进型”SEP事件。9月29日22点12分（UT)），部署在太阳-地球引力平衡点L1处的SOHO卫星观测到了全晕CME事件的爆发。10月2日2时（UT）左右，ACE 卫星观测到CME 激波到达地球附近。GOES15卫星所搭载EPEAD（Energetic Proton Electron and Alpha Detector）载荷的观测显示，在此次事件之前20多天的时间内，高能粒子通量处于稳定的低水平背景状态，在事件结束后的20多天内没有观测到典型的太阳高能粒子事件。这表明本文所研究的这一事件在时间上具有较好的孤立性，可以排除多个事件相互耦合所造成的复杂影响。SOHO卫星的观测显示，在此次全晕CME事件发生的前后一段时间内，没有其他全晕CME 的发生。

如前所述，模型设置的网格为1500×360，为简单起见，模型在计算背景太阳风参数时假设内边界上所有点处的内边界条件都是轴对称的，即同一物理量的数值在内边界所有的点上都是相同的。通过ACE卫星的观测获得1 AU 附近（即地球附近）的背景太阳风参数以及采用一维多方模型计算得到的内边界处的背景太阳风参数分别如表1所示。其中在计算磁场分量时，假设磁场在z方向的分量很小，因此磁场分量的计算结果简化为二维平面中的2个分量，即径向分量和切向分量。

表1 背景太阳风参数Table 1 Background solar wind parameters

本文利用冰激凌锥模型[31]来反演CME参数。首先使用人工方法对SOHO卫星日冕仪的白光日冕图像进行CME 前沿的信息提取。图1分别展示了对不同时刻的C3图像进行人工CME信息提取的过程，其中实心点为人工标记的CME 前沿。

图1 人工CME 前沿标记Fig.1 Artificial CME frontier marking

进而通过计算，由反演程序输出了图2所示的拟合曲线，并获得了表2中的CME 参数。图2显示了冰激凌锥模型的拟合曲线、CME 抛射速度和CME 抛射的角宽度。

图2 冰激凌锥模型拟合曲线Fig.2 Fitting result of the ice-cream cone model

表2 冰激凌锥模型反演结果Table 2 Inversion result of the ice-cream cone model

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

2 数值模拟结果

2.1 iPATH模型实例模拟结果

利用上述方法得到的模型输入参数，以地球位置为观测点进行了数值模拟计算。在模型中设定CME 中心方位角为90°，并通过冰激凌锥模型反演的源区位置计算了地球在坐标系中相应的位置。图3为归一化数密度的数值模拟结果，分别展示了4 个不同时刻的CME 和激波在黄道面内传播的过程中磁力线弯曲的状态和演化。图3 中：绿色标记点为地球位置；螺旋实线为磁力线，其中绿色实线为穿过地球的磁力线；颜色棒代表了归一化的密度值，颜色越深代表数值越大。由图可见，观测点与激波波前处磁力线连接的位置对太阳高能粒子的模拟结果有着非常重要的影响；激波在传播的过程中，逐渐由准平行激波转化为准垂直激波，同时连接点的位置也从激波边缘逐渐向激波中心移动。从图3(c)可以看出，激波在10月2日13：42：45时刻已经越过了地球所在位置，这与ACE 所观测到的激波到达时间基本相符。

图3 CME 和激波在黄道面中的传播Fig.3 The simulated propagation of the CME and the shock in the ecliptic plane

图4对比了此次太阳高能粒子事件的数值模拟结果与GOES15卫星的观测结果，其中：红色、蓝色、黑色的曲线分别代表能量≥10 MeV、≥50 MeV和≥100 MeV 高能质子的积分通量；实线和虚线分别代表GOES15的观测结果和数值模拟的结果。由于CME 爆发在太阳表面的附近，而模型的内边界在0.1 AU 处，所以模型的计算结果忽略了CME从太阳表面附近传播到0.1 AU 处的过程。在对比时，通过反演得到的CME 抛射速度估算了CME传播到0.1 AU 处的时间，从而对模拟结果中高能粒子事件的开始时间进行了简单的校正。由图可见：数值模拟结果中＞10 MeV 高能粒子通量与观测结果吻合较好，其变化趋势与观测结果较为相近；对于能量＞50 MeV 和＞100 MeV 高能质子的积分通量而言，模拟结果虽然能够大致体现出与观测结果相应的变化趋势，但其通量高于观测值，尤其是能量＞100 MeV 的高能质子积分通量差异显著。造成这种结果的可能原因之一是，由于在模拟中对粒子注入率的处理与实际物理过程有很大差异，导致较高能段的粒子注入率偏高。另外一种可能的因素是二维数值模拟不能很好地模拟三维空间中的CME传播，导致激波压缩比计算值偏高。

图4 GOES卫星观测的SEP 事件单向积分通量与模拟结果对比Fig.4 Comparison between the observed unidirectional integrated flux obtained by GOES satellite and the simulated results

2.2 不同CME 抛射速度条件下的集合模拟结果

限于观测手段，锥模型拟合的CME 的抛射速度和主抛射方向往往和实际情况相差较大。为了研究不同的CME抛射速度对iPATH模型数值模拟结果的影响，本文分别选取了CME 抛射速度分别为900、1000、1100和1300 km/s进行数值模拟，并将模拟结果与前述模拟结果进行对比。图5显示了同一时刻不同CME抛射速度条件下的CME/激波状态，从图5(a)到图5(d)，CME抛射速度依次增大，其中图5(c)为2.1节中的模拟结果，与图3(c)相同。可以看到，在同一时刻CME抛射速度越高，激波波前距离太阳越远。图6显示了不同CME抛射速度条件下能量＞10 MeV 高能质子的单向积分通量。可以看到，积分通量曲线在不同CME抛射速度下有较为明显的差异：在事件开始的前段，积分通量呈现出随着CME 抛射速度的增大而增高的趋势；在事件后段，积分通量的差异开始收缩。

图5 不同CME 抛射速度条件下的CME/激波对比Fig.5 Comparison of CME/shock propagations of different initial CME speeds

图6 不同CME 抛射速度条件下＞10 MeV 的高能质子通量模拟结果对比Fig.6 Comparison of simulated results of energetic flux of particles of more than 10 MeV with different initial CME speeds

进一步考查不同CME 抛射速度条件下高能粒子能谱的特征。图7中所示为不同CME抛射速度条件下整个事件的高能粒子能谱分布以及与之相对应的分段线性拟合的结果，其中：不同颜色的实心点代表不同抛射速度，直线代表对能谱进行分段拟合的结果。可以看出，CME抛射速度越大的高能粒子能谱越高。

图7 不同CME 抛射速度条件下高能质子全事件积分能谱分布及拟合结果Fig.7 Event-integrated power spectra and fitting results for different initial ejection speeds

为了定量显示高能粒子能谱间的差异，使用了形如y=kx+b的函数对高能粒子能谱进行了最小二乘拟合，并给出了与之相应的分段拟合结果，如表3所示。表中给出了拟合直线的k和b，其中下角标1、2分别表示第1段和第2段的拟合直线。结果显示粒子能谱间具有较大差异。

表3 不同CME抛射速度条件粒子能谱拟合结果Table 3 Fitting result of the particle spectrum for different CME ejection speeds

2.3 不同内边界背景太阳风温度条件下的集合模拟结果

一般情况下，内日冕的太阳风温度无法直接测量。为进一步研究不同内边界背景太阳风温度对iPATH 模型数值模拟结果的影响，采取与2.2节相似的研究方法。在此次事件中，我们通过一维多方模型估算的内边界背景太阳风温度为0.4287 MK。本文分别选取此数值的0.5倍、0.75倍、1.5倍和2 倍作为内边界温度条件进行数值模拟，并将结果与前述模拟结果进行对比。图8所示为同一时刻不同内边界背景太阳风温度条件下的CME 激波状态，由图8(a)到图8(d)温度依次增大，其中图8(c)为前述实例模拟的结果。可以看到，在同一时刻背景太阳风温度越高，激波波前面距离太阳越远，但这种变化相对于图5来说明显更小。这说明背景太阳风温度的增高会引起CME 速度的微弱变化。图9 所示为不同内边界背景太阳风温度条件下能量＞10 MeV高能质子的单向积分通量。可以看到：积分通量曲线在不同内边界背景太阳风温度下的差异很小；在事件开始的前段，积分通量呈现出随着内边界背景太阳风温度增大而有微弱增高的趋势，在事件后段则呈现出与前段相反的趋势。

进一步考查不同内边界背景太阳风温度条件下高能粒子能谱的特征。图10所示为不同内边界背景太阳风温度条件下的高能粒子能谱分布以及与之相对应的分段线性拟合的结果，其中：不同颜色的实心点代表内边界处不同的背景太阳风温度，直线代表对能谱进行分段拟合的结果。可以看出，内边界背景太阳风温度的变化对高能粒子能谱的影响很小。为了定量地显示高能粒子能谱间的差异，本文仍然使用形如y=kx+b的函数对高能粒子能谱进行了最小二乘拟合，并给出了与之相应的分段拟合结果，如表4所示。结果清晰地表明内边界背景太阳风温度的变化对SEP事件中高能粒子能谱的影响很小。

图8 内边界不同背景太阳风温度条件下的CME/激波对比Fig.8 Comparison of CME/shock propagation at different background solar wind temperatures

图9 内边界不同背景太阳风温度条件下＞10 MeV 的高能质子通量模拟结果对比Fig.9 Comparison of the simulated results of energetic flux of particle of more than 10 MeV at different background solar wind temperatures

图10 内边界不同背景太阳风温度条件下的高能质子全事件积分能谱分布及拟合结果Fig.10 Event-integrated power spectra and fitting results for different background solar wind temperatures

表4 不同背景太阳风温度条件下粒子能谱拟合结果Table 4 Fitting result of the particle spectrum at different background solar wind temperatures

3 结束语

本文使用iPATH 模型对2013年9月30日的一次孤立SEP事件进行了数值模拟。数值模拟结果和GOES卫星观测结果的对比显示：数值模拟的能量＞10 MeV 的高能粒子的通量和观测值吻合较好，而能量＞100 MeV 的高能粒子通量的模拟结果高于观测值。造成这种情况的可能原因之一是不同能量的高能粒子在加速过程中经历不一样的粒子注入过程，或者说注入率不一样。另外一种可能的原因是我们在计算过程中采用的是二维模型，因此对激波参数的计算和实际的三维情况存在偏差。另一方面，基于SOHO卫星白光日冕观测图像，通过冰激凌锥模型反演CME参数，然而模型反演的CME 参数往往和实际情况相差较大。鉴于SEP事件的形成过程受到CME 初始抛射速度的影响较大，为了定量地研究CME 初始抛射速度带来的影响，进一步开展了不同CME 初始抛射速度下的集合模拟试验。集合模拟的结果表明更高的CME 初始抛射速度会引起更强的SEP事件。

我们的预报模型采用了一个较为简单的一维多方模型来估计内边界处背景太阳风的参数。一般而言，内边界处背景太阳风的速度、密度和磁场的估算较为容易，而温度的估算则相对困难，一维多方模型所估算的温度往往高于实际的温度。此外，以往的文献基本上没有关于背景太阳风温度对于SEP事件影响的相关讨论。因此，我们开展了不同背景太阳风温度条件下的集合模拟研究。数值计算的结果表明：内边界背景太阳风温度的变化对于地球附近SEP事件中高能粒子通量和能谱的影响几乎可以忽略不计。该结论在空间环境态势感知中具有重要的应用意义。在为各类空间探测计划提供保障服务、实际进行太阳质子事件预报的过程中，即使我们不能非常准确地获取靠近太阳的内边界处背景太阳风的温度，也不会显著地影响预报的结果。

致谢

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