2 AU以内的“渐进型”太阳高能粒子事件模拟

敖先志，刘四清,2，沈 华，王晶晶，胡骏翔，李 刚

（1.中国科学院国家空间科学中心，北京 100190；2.中国科学院大学，北京 100190；3.美国阿拉巴马大学汉茨维尔分校空间科学系，汉茨维尔 35899）

引 言

太阳高能粒子事件是影响近地空间（1 AU）以及深空辐射环境重要的甚至是主要的因素之一，对航天器、人造地球卫星，以及宇航员的安全等构成巨大的威胁。大规模的太阳高能粒子事件会造成严重的空间灾害天气。因此，太阳高能粒子事件的观测、产生形成机制以及相应的预报模式是当前国际空间物理学界的最前沿课题之一。

太阳高能粒子事件往往可以分为两种类型：“脉冲”（impulsive）型和“渐进”（gradual）型[1]。“脉冲”型事件中的高能粒子表现在观测上的特征是仪器测量到的高能粒子通量有个突然的、急剧的增强，然后又快速下降，持续时间常常在一个或者数个小时左右，就好像一个单脉冲一样；“渐进型”事件中的高能粒子其典型观测特征是升高后的粒子通量缓慢、渐进地下降，事件持续的时间很长，可以达到2 天以上。这两种形态的事件有时候会混合在一起[2-4]。在实际观测中，“渐进”型的太阳高能粒子事件占大多数[5]。

国外对于太阳高能粒子事件的观测和理论研究开展得比较早[1,6]。一般认为，“渐进型”事件中的高能粒子（主要成分为质子）是由日冕物质抛射驱动（CME-driven）的激波加速产生[7-11]。这种粒子加速机制通常被称为激波扩散加速（Diffusive Shock Accelera‐tion，DSA）机制，其本质是一阶费米加速。这种加速机制可以轻而易举地把粒子加速到＞100 MeV/nucleon。2000年，Zank 等提出了基于DSA 加速机制的太阳高能粒子加速和传播模型：PATH（Particle Acceleration and Transport in the Heliosphere）模型。这个模型在过去的10多年里得到了不断的发展和完善[12-20]。胡骏翔等进一步扩展了前人的研究[21-22]，建立了iPATH（improved PATH）模型，将激波的演化从一维扩展到了二维，从而能够考虑方位角和垂直扩散系数对粒子加速和传播的影响。

本文利用iPATH模型对发生于2014年04月18日的CME 所引起的太阳高能粒子事件实例进行了数值模拟，并将地球附近（1 AU）的数值模拟结果和卫星观测进行了对比；在此基础上加以扩展，进一步给出了深空环境中不同位置处太阳高能粒子通量的模型计算结果。

1 iPATH模型简述

iPATH模型是在PATH模型的基础上进行的扩展。iPATH模型的运行主要包括3个步骤：①生成背景太阳风环境；②模拟CME 驱动的激波并计算粒子在激波附近的加速过程；③计算被激波加速的高能粒子在行星际空间中的传播。

iPATH对背景太阳风的模拟使用了开源数值模拟程序包ZEUS3D[23]来求解磁流体力学方程组为

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

尽管ZEUS3D 具有求解三维磁流体力学（Mag‐netoHydro Dynamics，MHD）方程组的能力，但是iPATH在当前只具有模拟二维条件下的粒子加速和传播能力，因而在利用ZEUS3D 程序包进行计算时，iPATH模型将其中一个维度的大小设为1。具体来说，ZEUS3D 程序包的设置采用了球极坐标（r，θ，φ）参考系，r，θ，φ分别为径向距离、仰角和方位角，太阳位于中心，θ=90°的平面是黄道面，iPATH模型的θ角取值被限定于只取90°一个值，从而使得三维模拟程序退化为二维数值模拟。这种方法的不足在于二维数值模拟没有三维数值模拟更准确，但是对于太阳高能粒子事件模拟来说可以减少计算资源、大幅提高计算速度，尤其是对于空间天气预警和预报而言，计算时间的减少，带来的是预警提前量的提升。

iPATH模型通过在内边界处引入持续一段时间的扰动来模拟日冕物质抛射。内边界处的扰动参数包括CME 速度、密度、温度、磁场以及扰动持续时间，CME的角宽度也是参数之一，用于描述CME的扰动速度分布的范围。为简单起见，iPATH 模型假设CME的速度分布为高斯分布。

对于发生于2014年4月18日的CME事件，本文所选择的初始扰动参数，尽量使得模拟结果符合日冕仪和1 AU 处的太阳风的实际观测数据。与背景太阳风数值模拟类似，iPATH 模型采用MHD 方程组来描述扰动在太阳风中的传播。随着CME 扰动的引入，CME驱动的激波开始在模型中传播和演化。ZEUS3D程序包采用了交错网格和迎风格式。模型中的初始磁场为简单的Parker螺旋场。

iPATH模型在计算激波参数时同步计算粒子的加速过程。对于粒子在激波附近经历加速过程的理论和数值模型，文献[12]～[22]有非常详尽的描述，本文不再重复叙述。文献[12]～[20]所发展的PATH模型在后期已经在一定程度上考虑了垂直扩散系数，但是鉴于PATH模型在处理激波时，采用的是一维模型，因此PATH 无法解决扩散系数在φ方向上的变化。iPA‐TH 模型从一维扩展为二维空间，扩散系数进而可以表现为φ的函数，或者更为明确地说，激波波前不同位置具有不同的扩散系数。

高能粒子的注入率对被加速粒子能谱的形成具有重要作用。注入率受到局地太阳风参数、湍流强度、激波参数、粒子能量动量等各种因素影响，也是目前高能粒子加速机制研究的一个重要对象。现阶段，在iPATH模型里粒子的注入率是高能粒子能量动量的一个函数。此外，iPATH模型假设：激波对粒子有加速作用，而粒子对于激波的反作用可以忽略。iPATH模型利用磁流体力学方程组来描述和模拟激波传播，然而磁流体力学是没有办法来刻画高能粒子加速过程对激波反作用的动力学特征的。这一点，有待于未来的研究工作来解决。

被加速的高能粒子从激波附近逃逸并沿行星际空间磁力线传播，传播过程中被传播路径上的湍流所散射。粒子的传播过程服从聚焦传输（focused trans‐port）方程，在iPATH 模型中利用蒙特卡洛（Monte Carlo）方法来求解粒子的传播过程[21]。

2 数值模拟结果

数个 GOES（Geostationary Operational Environ‐mental Satellites）卫星从 2014年04月18日—2014年04月22日均观测到了强烈的太阳高能质子增强事件。这次太阳高能粒子事件来源于一次大的日冕物质抛射事件所驱动的激波。在2014年04月18日，太阳活动区AR12036 爆发了一次M7.3 级的耀斑，X 射线的流量在当天13：03UT左右达到了最大值。耀斑爆发大约20 min 后，位于日地引力平衡点L1 处的太阳和日球层探测器（Solar and Heliospheric Observatory，SO‐HO）卫星所搭载的C2 和C3日冕仪观测到了高速全晕CME。先进的成分探测器（Advance Composition Explorer，ACE）卫星在 20日大约 10：20UT 左右观测到了这一次CME爆发所驱动的激波。这次全晕CME事件所引发的太阳高能粒子事件是一起孤立事件，事件前后高能粒子的通量均无由于其它CME 事件所带来的明显扰动。GOES-15 卫星的观测表明在事件之前7 天和之后9 天的时间内，能量大于2.5 MeV 的质子通量保持在背景大小。因此，这次事件非常适合进行数值模拟研究。实际上，在此次CME 爆发的同一天大约07：24UT 左右爆发了一个小的CME，其速度和角宽度经估算分别为大约387 km/s 和84°。同时期ACE 卫星观测到的太阳风速度逐渐从大约400 km/s增加到700 km/s，因而这一次小的CME 事件不大可能会引起SEP事件。另一方面，STEREO的2颗卫星也观测到了这次小的CME 事件。STEREO-A 与地球的夹角为165°，STEREO-B与地球的夹角为156°，这两颗卫星上搭载的日冕仪观测图像结合SOHO卫星的日冕仪图像，我们推断07：24UT 左右爆发的小型CME 事件的抛射方向远离黄道面，与日地连线接近垂直。故而，这次事件不会对本文的数值模拟造成非常显著的干扰。

此次事件数值模拟计算域的物理尺度为

其中：r的单位为AU；φ的单位为角度。

网格大小设置为黄道面上的1 500×360 个格点。在r方向和φ方向上均为均匀网格划分，因此计算域的内边界是一个半径为0.1 AU 的圆。在计算背景太阳风参数这一步中，为简化起见，采用了轴对称设置，即内边界所有网格点上同一个物理量的数值是相同。数值模拟所使用的内边界背景太阳风参数见表1所示。图1给出了背景太阳风的模拟结果。图中横坐标是离太阳的距离（以AU 为单位），（a）是归一化的太阳风数密度，即质子数密度乘以距离的平方；（b）是太阳风的速度。在1 AU 处，数值模拟的太阳风数密度大约为5.9/cm3，速度大约为468 km/s。同时期的ACE卫星在1 AU的观测数据表明，太阳风数密度围绕6/cm3上下波动，速度围绕500 km/s上下波动。因此，采用表1中的内边界条件数值模拟得到的结果是符合实际太阳风状态的。

表1 数值模拟内边界参数Table1 Inner boundary conditions

图1 背景太阳风参数数值模拟结果。上：归一化的太阳风质子数密度；下：太阳风速度。横坐标是r方向上的距离Fig.1 Numerical simulation of the background solar wind conditions

中国科学院国家空间科学中心所属空间环境预报与研究中心（Space Environment Prediction Center，SEPC）部署了一套日冕物质抛射自动识别软件程序。这套软件利用Hough 变换和J-maps[24-25]方法从SOHO卫星日冕仪的白光日冕图像中自动检测和识别CME，并进一步采用冰激凌锥模型[26]来反演CME 的速度和角宽度等参数。在某些情况下，例如全晕CME，自动识别和反演程序可能会把1个CME识别成2个或者多个CME，SEPC 的CME 识别软件同时提供人工干预接口以除去明显不合理的结果。

图2是计算机程序自动识别和人工干预识别对2014年04月18日的CME爆发进行识别的对比图。图2（a）是自动识别的结果，图2（b）是人工干预识别的结果。自动识别将此次全晕CME 识别为2 个CME，本文在识别CME时进行了人工修正。将CME识别结果利用冰激凌锥模型进行反演，结果见图3。图3（a）是将图2（b）的结果输入冰激凌锥模型反演的结果，图3（b）是全自动反演的结果。显然，图3（b）结果不是真实的。因此，本文在后续的数值模拟中采用图3（a）的结果：速度1 290 km/s，角宽度120°。此外，数值模拟所采用的CME 持续时间为1.5 h，其它的CME参数见表1。

图2 SOHO卫星白光日冕图像CME识别与提取Fig.2 CME detection and recognition based on SOHO coronagraph images

图3 冰激凌锥模型反演CME参数Fig.3 CME parameters derived by Cone model fitting method

图4所示为CME/激波在黄道面内传播的数值模拟结果。每个小图中，太阳位于中心，以白色表示；绿色代表地球所在的位置，蓝色、红色和黄色代表另外3个假想的观测点，其中红色和黄色观测点的方位角分别为30°和5°，地球和蓝色观测点的方位角为90°；黑色曲线代表磁力线；颜色棒表示的是速度的大小，单位为km/s。观测点和激波波前的磁力线连接位置（Connecting with the OBserver Point，COB‐Point）在太阳高能粒子事件中起着非常重要的作用。这几个观测点在研究磁力线联通与SEP事件之间的关系具有相当的代表性。例如，红色观测点在爆发开始时，磁力线连接在激波最强的中心区域，随着时间的推移，连接点（COBPoint）逐渐向激波右侧边缘移动；黄色观测点的磁力线连接点（COBPoint）起初位于激波左侧边缘，随激波传播逐渐向激波中心区域移动；蓝色观测点在开始时刻几乎不与激波波前具有磁力线连接。图3中径向距离的单位为AU，黑色实线圆为1 AU 的地球轨道。图4中可以直观地看到，激波法向和磁力线的夹角随激波波前位置不同而变化。图4中每一个小图，激波中间部分压缩比大，向两侧，压缩比逐渐减小。因此，不同位置处被加速的粒子也具有不一样的能谱。

高能粒子的注入率在SEP事件数值模拟中扮演着非常重要的作用，SEP 事件中的粒子通量均与之相关。假设粒子的注入率为[20]

图4 CME/激波在黄道面内传播Fig.4 The simulated propagation of the CME/shock in the ecliptic plane

其中：θ为激波法向和磁力线的夹角；E0是质子在平行激波条件下的注入能；是质子在斜激波条件下的注入能，它是θ的函数，( -δ)是种子粒子的能谱分布函数的幂指数。

在深空环境下，很少有机会能够直接将观测结果和模拟结果进行对比。本文利用地球轨道的观测结果来确定数值模型的参数，再以同样的参数来模拟深空环境下的SEP事件。如前所述，假设粒子注入率中的χ= 2%，数值模拟所得的2014年04月18日SEP事件在1 AU地球附近的高能质子的单向积分通量见图5。

图6是数值模拟所得到的地球附近的高能质子的单向微分通量，图6给出了6个不同能段粒子的单向微分通量。总而言之，对于2014年04月18日的SEP事件在地球附近的表现所进行的数值模拟基本上能反映实际的观测结果，本文利用地球观测的结果所确定的数值模拟参数，将模拟范围扩大到1 AU 以外，进一步对深空环境下的SEP事件进行模拟。

图5中虚线是地球同步轨道GOES卫星的观测结果，实线是iPATH模型的模拟结果。红色、蓝色和绿色曲线分别代表能量大于10、50、100 MeV 的高能质子的单向积分通量。以大于10 MeV的高能质子为例，数值模拟的结果在事件发生后的前12 h内和观测较为吻合，之后我们模拟的积分通量开始下降，与之相比观测值仍然维持了很长时间才开始下降。通过与GOES卫星观测的高能质子的单向微分通量以及软X射线的观测进行对比发现，在事件发生的大约前12～24 h，数值模拟的结果相对于观测结果偏小，其原因可能与行星际太阳风结构有关。在大约27 h后，太阳有一次小的耀斑爆发，耀斑爆发产生的入射粒子可能作为种子粒子注入激波附近进一步增强了这次SEP事件中＞10 MeV 能量粒子的通量，但是其对更高能的粒子＞50 MeV 的通量影响不大。事件发生后的前12 h内，对于＞50 MeV和＞100 MeV的高能粒子，本文数值模拟的结果偏大。

图5 地球附近SEP事件单向积分通量观测和模拟对比虚线为同步轨道GOES卫星观测值，实线为数值模拟结果Fig.5 Comparison between the observed unidirectional integral flux（dashed lines）by GOES satellite and the simulated results（solid lines）

图6 地球附近SEP事件单向微分通量的数值模拟结果Fig.6 The simulated unidirectional differential flux at near-Earth space

图7所示为1 AU距离30°处的观测点高能质子单向积分通量的数值模拟结果。由于该位置处的磁力线在CME 爆发开始时刻就连接于激波波前中心区域，因而高能粒子的通量迅速达到最高值，随着激波快速向外传播，磁力线在激波波前的连接点（COBPoint）也迅速地移动到激波波前边缘，在事件发生后大约15 h后，各个能段高能粒子的单向微分通量渐渐回归到背景大小（图8）。

图7 红色观测点位置处的SEP事件单向积分通量数值模拟结果Fig.7 The simulated unidirectional integral flux at location of the red spot

图8 红色观测点位置处的SEP事件单向微分通量数值模拟结果Fig.8 The simulated unidirectional differential flux at location of the red spot

图9～10是位于黄道面上90°，1.5 AU 处高能粒子的数值模拟结果。由于事件开始时刻观测点的磁力线并未很好地连接在激波波前，因此起初在该位置较低能段高能粒子的通量没有抬升，而较高能段的高能粒子由于横向的扩散，其通量仅有较小的抬升，在最初的5 h 左右，＞50 MeV 的高能质子通量的数量级从10-2上升为10-1。随着激波向外传播，该处的磁力线在激波波前的连接点（COBPoint）向激波中间区域移动，低能段的高能粒子通量也开始抬升，但是由于激波在传播过程中逐渐减弱，对高能粒子的加速作用也降低，因而高能段的高能粒子通量的抬升也不大，从图10中可以看到，18 MeV的高能粒子单向微分通量最高到了大约0.4 个/（cm2s sr MeV），而1.1 MeV的高能粒子的单向微分通量最高则上升了到了接近100个/（cm2s sr MeV）。

图9 蓝色观测点位置处的SEP事件单向积分通量数值模拟结果Fig.9 The simulated unidirectional integral flux at location of the blue spot

图10 蓝色观测点位置处的SEP事件单向微分通量数值模拟结果Fig.10 The simulated unidirectional differential flux at location of the blue spot

图11～12是位于黄道面上 5°，1.5 AU 处的 SEP数值模拟结果。在事件初始时刻，此处磁力线在激波面上的连接点（COBPoint）靠近左侧，伴随连接点逐渐向激波中心区域移动，高能粒子的通量开始抬升，到达最高点后开始下降。随着连接点逐渐向激波右侧边缘移动，激波的加速作用很快降低，在大约25 h 后观测点处的高能粒子通量基本恢复到背景水平。

图11 黄色观测点位置处的SEP事件单向积分通量数值模拟结果Fig.11 The simulated unidirectional integral flux at location of the yellow spot

图12 黄色观测点位置处的SEP事件单向微分通量数值模拟结果Fig.12 The simulated unidirectional differential flux at location of the yellow spot

3 结 论

本文利用iPATH模型对发生于2014年04月18日爆发的CME 所驱动的激波引起的行星际空间太阳高能粒子事件进行了数值模拟。数值模拟局限于2个天文单位距离以内的二维黄道面内。本文在模拟中选取了4个具有代表性的观察点，其中一个观察点是地球所在位置。通过对比模拟结果和GOES卫星的观测结果来确定合适的数值模拟参数，并以此为基础扩展到了更遥远的深空。结果表明：黄道面内不同距离不同角度的位置，在一次SEP事件中可以经历完全不同的过程。

现阶段太阳高能粒子事件的数值模拟还较为粗略，由于整个SEP事件的物理过程非常复杂，iPATH模型做了必要的简化，所采用的方法和结果不但揭示了深空环境中SEP的物理过程，而且也为深空环境的高能粒子辐射分析提供了参考。深空探测计划必然需要考虑项目实施过程中来自太阳高能粒子事件的影响，本研究工作所获得的成果可以在这方面得到应用。