0 引言

近地空间的太阳风参数预报具有重要的科学研究意义和实际应用价值(Feynman and Gabriel, 2000).三维MHD数值模拟是太阳风参数预报的重要手段，它能够提供太阳风在日地空间的分布和演化，提前给出太阳风参数在近地空间的变化.近二十年来，太阳风的MHD数值模拟取得了很大的进展，国内外已经发展了诸多数值模型(例如Riley et al., 2001; Odstrcil, 2003; Detman et al., 2006; Feng et al., 2007, 2010; Tóth et al., 2012; Hayashi, 2012; Wiengarten et al., 2014; Shiota et al., 2014等)，这些模型已经从定性的研究逐步过渡到集成预报模式的应用以及事件的定量分析研究，关于这方面的详细进展可参考冯学尚等(2011)、Wu和Dryer(2015).

太阳风的MHD数值模型可以分为两类：一类是基于物理过程的模型，以日冕底部为内边界，在方程中添加太阳风的加速和加热等物理过程的数学描述，以产生更切近实际的太阳风背景，例如CSEM模型(Tóth et al., 2012)、SIP-CESE-MHD模型(Feng et al., 2007)等；另一类是基于经验模型，这类模型一般将内边界设置在阿尔芬临界点之外，忽略日冕的具体物理过程，用经验模式来确定边界条件，模拟太阳风在行星际空间的演化过程，例如WSA/ENLIL模型(Odstrcil, 2003)等.由于太阳表面β≪1(β为等离子体热压与磁压之比)，前者在日冕部分消耗了相当多的计算资源；后者以经验关系作为边界条件，省去了日冕的计算资源的开销，计算效率得到大幅提高.Owens等(2008)详细对比了两类模型，结果显示基于经验模式的数值模型在再现太阳风的大尺度结构特征方面的效果并不亚于基于物理过程的数值模型.因此，无论是从时效性还是准确性方面考虑，基于经验模式的数值模型在太阳风预报方面都有重要的应用价值.

基于经验模式的数值模型能否取得与观测符合更好的结果很大程度上取决于底部所采用的边界条件.Owens等(2008)指出通过调试边界条件，可以很大程度上改善模拟结果.目前已发展的基于经验模式的数值模型的边界条件的取法各有特点，关于这些模型的具体做法可以参考张嫚和周玉芬(2014).现有的模型普遍存在的问题是给定的边界条件往往不能兼顾所有的物理量，直接导致模拟结果不能反映行星际空间的真实状态，所以有必要改进边界条件的确定方法及其自由参数的取值.为了能反映太阳风参数随太阳活动周的变化，本文基于冯学尚开发的TVD-COIN模型，提出了一组完备自洽的边界条件，建立了行星际太阳风的三维MHD数值模型，模拟可以呈现不同太阳活动相期间的太阳风结构.背景太阳风在一定程度上影响着地磁活动(Richardson et al., 2002)，并且提供更准确的背景太阳风的状态对于提高驱动大型地磁暴的行星际日冕物质抛射到达地球时间的预报非常重要(Wu et al., 2005)，因此本文重点考虑用建立的模型模拟背景太阳风.第1节介绍三维MHD数值模型，第2节结合WSA模型以及太阳风卫星观测特征提出一组新的边界条件，第3节利用上述模型分别模拟2007年和2016年的背景太阳风，并将模拟的结果和近地卫星ACE/WIND的观测数据对比.

1 数值模型

COIN-TVD模型是一套组合数值格式的三维太阳风模型(冯学尚等, 2002；Feng et al., 2003, 2005；Shen et al., 2007, 2009).原模型分为日冕和行星际两个求解区域：在日冕区域先用TVD Lax-Friedrich差分格式计算出新时间步的值, 再将磁场部分用MacCormack型格式的计算结果替换掉；在行星际空间用MacCormack格式求解与时间无关的定态MHD方程组，这样就完成了一步求解.Shen等(2007, 2009)对模型进行了修改，在全部区域采用TVD Lax-Friedrich格式求解三维时变MHD方程组，从而可以高效稳定地实现事件模拟.利用该模型模拟具体事件，结果再现了CME传播过程中在经度方向上的偏转, 1AU处的总磁场、密度、温度和速度剖面与观测剖面大体一致，但是部分物理量的变化不够理想，因此模型仍然有需要改进的地方(冯学尚等, 2011).针对模拟行星际太阳风三维时变结构的需要，本文进一步完善了模型，接下来将对改进后的模型进行简单介绍.

模型求解的是三维理想MHD方程组，物理量包括等离子体密度ρ，压强p，速度矢量V和磁场矢量B.由于太阳风和行星际磁场在行星际随着太阳的自转而刚性共转，所以选取共转坐标系可以方便底部边界条件的处理.在共转坐标系下，守恒形式的理想MHD方程组为：

(1)

(2)

(3)

(4)

其中G为万有引力常数，M_s为太阳质量，μ₀为真空磁导率.动量方程中f = -ω×[ω × r + 2ω × V]是因坐标系选取而引入的离心力，ω为太阳自转的角速度.忽略太阳表面的较差自转(Schröter, 1985)，在模型中使用固定的自转周期T_CR=25.38天，相应角速率ω=14.38°/天.考虑到数值格式的稳定性，在方程组中使用了压力方程替代能量方程(Odstrcil, 2003).压力方程中γ为多方指数，它的大小与太阳风的速度结构无关，取值为1.46(Totten et al., 1995).

在数值计算过程中，磁场的散度一般不能保证精确为零.如果不及时加以控制，磁场散度就会迅速增大，进而引起沿着磁力线方向上的非物理流动，甚至导致程序崩溃.为了消除磁场的散度，模型中采用扩散法(Rempel et al., 2009; Van Der Holst et al., 2007; Feng et al., 2010).这种方法相当于在磁感应方程中添加扩散项，使数值计算产生的磁场散度在计算区域内迅速扩散开来，从而避免了其在某一区域累积.扩散法的迭代形式为

(5)

其中，系数μ反映扩散的快慢，当使用中心差分离散时，根据稳定分析得到μ∈(0, 2)；Δx与网格大小有关，在球坐标系下(Δx)² = .数值模拟表明，在每一时间步之后，使用不超过10次的迭代可以满足磁场散度的相对误差.

模型中计算区域的底部边界R_b=21.5R_s(太阳半径)，即0.1AU处.在这一位置处，无论是快速流，还是慢速流，都基本完成加速过程，一般情况下满足超声速以及阿尔芬速的条件，因此所有的特征波都向外传播，计算区域内的扰动不会通过边界传回日冕.此外，高速流与低速流之间的相互作用不是很强，等离子体携带的磁场在日冕经过充分的膨胀后也大体上沿着径向方向(Zhao and Hoeksema, 2010)，所以这一位置比较方便边界条件的给定.在径向方向上采用等比网格，求解范围包括0.1AU到1.1AU，网格大小从0.27R_s逐渐增加到2.15R_s, 网格总数为200个.根据边界条件的特征分析方法，在内边界上所有特征波指向计算区域内，需要给定全部物理量；在外边界上所有特征波指向计算区域外，各物理量根据内点线性外推得到.在球面上采用六片网格系统(Feng et al., 2010)，它由六个完全相同低纬区域构成，不同片之间相互交叠，覆盖整个球壳区域，每片区域内部网格的物理量通过数值格式求解，边界上的网格的物理量则通过与其重叠的相邻片上的内部网格插值得到，网格的空间分辨率为1°×1°.六片网格系统的优点是它避免了传统球坐标系中极区的奇异性和收敛性问题(Wu et al., 2015)，并且大幅增加时间步长，提高了计算效率.

2 边界条件

如前所述计算区域的内边界位于超声速区域，因此在内边界上的物理量需要全部给定.由于目前缺乏直接的观测，一般的做法是借助经验模型获得太阳风的速度和磁场的径向分量的，然后结合其他假设导出密度和温度以及速度和磁场的切向分量，最终构成一组完备的边界条件.这一过程的出发点是太阳光球层的观测磁图，目前有多种地基或者卫星观测磁图可供选择，本文采用的是GONG的磁图(http://gong.nso.edu/).结合太阳风的卫星观察特征，本文提出了一套新的边界条件的确定方案.

在现有的模型中，内边界上磁场的径向分量B_r大多通过势场源表面+Schatten电流片(PFSS+SCS)模型反演观测的太阳磁图得到(Schatten et al., 1969; Altschuler and Newkirk, 1969; Schatten, 1971). PFSS+SCS模型不仅能很好预报行星际磁场的极性(Arge et al., 2003)，而且利用它得到的行星际磁场径向分量也和Ulysses观测到的与纬度无关的特征相符(Smith and Balogh, 1995).PFSS+SCS模型可以利用球谐函数展开法解析求解，本文设置模型外边界在0.1AU处，即太阳风数值模型的内边界所在的位置，勒让德多项式最大的级数L_max为11.考虑到现有的MHD模型普遍低估行星际磁场的通量(Stevens et al., 2012; Gressl et al., 2014)，本文采用1AU的观测来限制磁场大小：

(6)

其中mean(B^1AU)是ACE/WIND观测的磁场的强度在过去的三个卡林顿周的平均值，在计算之前根据Richardson and Cane的列表(http://www.srl.caltech.edu/ACE/ASC/DATA/level3/icmetable2.htm)剔除ICME(Interplanetary Coronal Mass Ejection, 行星际日冕物质抛射)的数据.由于在1AU处mean|B_φ|≈mean|B_r|, 所以再除以.令B₀= (B^1AU)，计算出B₀后再根据磁通量守恒并结合PFSS+SCS模型即可得到内边界上B_r的分布.

太阳风的径向速度是模拟中最重要的参数，这是因为在行星际背景太阳风主要的动力学过程是高速流和低速流之间的相互作用，速度的分布将很大程度决定其他所有太阳风参数的演化.在本文中，模型的内边界上太阳风的径向速度由WSA经验模型得到(Arge et al., 2003).Wang和Sheeley(1990)发现L1点的太阳风的速度和日冕磁通量管的膨胀因子f_s反相关；Riley等(2001)行星际等离子体数值模拟中，将内边界上的太阳风的速度和相应的磁力线在光球层的足点距离冕洞边界的最小角距离θ_b联系起来.Arge等(2003)结合上述两项工作的成果提出了WSA太阳风速度的经验关系，在21.5R_s太阳风的径向速度可以写成以下一般形式：

(7)

其中V_s，V_f分别限定太阳风最低速度和最高速度，其大小可以根据观测调节；f_s和θ_b可以根据PFSS+SCS模型计算得到.除此之外还有4个自由参数，a₁体现的是f_s的影响程度，a₂至a₄则反映θ_b的影响程度.WSA经验关系中的自由参数不是一成不变的，它们的最佳的取值随着时间变化(Riley et al., 2015).本文设定a₃和a₄均为1，仅保留a₁和a₂两个自由参数.考虑到模型在行星际空间对太阳风有小幅度的加速，利用上述关系式得到速度分布后需再减去50 km·s^-1(McGregor et al., 2011).

在行星际空间高速流与低速流的相互作用的过程中动压占据主导地位，给定合理的密度分布也是十分重要的.观测表明太阳风的密度并非独立变化，它与太阳风速度呈反相关，这意味着内边界上的密度不能任意给定.关于密度和速度二者之间的定量关系，不同观测数据得到的结果存在区别(例如Burlaga and Ogilvie, 1970; Steinitz and Eyni, 1980; Phillip et al., 1995).最近，Chat等(2012)根据Helios, Ulysses和Wind的数据发现太阳风的能量通量(在1AU主要包括太阳风等离子体的动能和重力势能)与速度、纬度无关，据此得到了密度的半经验关系式.本文根据Chat等的分析结果，在模型中将太阳风的能量通量设为常数，数密度的具体计算公式为：

(8)

其中N₀为1AU处速度V₀=750 km·s^-1的太阳风的数密度，其大小根据ACE/WIND观测的太阳风在过去的三个卡林顿周的平均能量通量(剔除ICME的数据后)确定.在其他MHD模型中，通常假设N∝V_r^-q，其中q取值为1或2.Detman等(2006)认为q应该介于1和2之间.实际上式(8)近似等效于q=5/3，它不仅满足黄道面内的观测，也与Ulysses的极区观测结果相符.

太阳风的质子温度与速度有着很强的正相关关系，因此可以采用类似计算密度的方法得到温度的分布.许多研究工作曾尝试拟合质子温度和太阳风速度在1AU处的定量关系(例如Burlaga and Ogilvie, 1973; Lopez and Freeman, 1986; Elliott et al., 2005; Chat et al., 2012).本文采用, 其中速度的单位为km·s^-1, 温度的单位为K.根据太阳风的Parker一维球对称多方模型(Parker, 1958)，假设太阳风速度在行星际沿着径向不变，那么温度将按照r^-2(γ-1)衰减.考虑到γ取值1.46，得到内边界上温度分布：

(9)

进一步假设质子和电子的温度相同T=T_p=T_e，可以得到内边界上的压强p=2NKT.

假设在静止坐标系下内边界上太阳风速度沿着径向方向，那么在经度和纬度方向上的速度分量可以设为零.通过坐标系变换，得到共转坐标系下的表达式：

(10)

由于太阳风的高电导率，行星际磁场冻结在太阳风等离子体中.在共转坐标系下，对于稳态结构，行星际磁场和太阳风速度的方向相互平行，据此可以计算内边界上磁场的其余分量：

(11)

至此，控制方程中所有物理量的边界条件都已经给出计算公式.为了尽可能使各个物理量的变化特征贴近观测，边界条件中保留了6个可调参数，其中径向速度包含4个，径向磁场和密度分别含1个，这些参数的取值范围可以参考表 1.特别地，表 2给出了2007年至2016年部分参数的具体取值.实际上，这些可调参数相当长的时期内(数个卡林顿周期)保持不变.在共转坐标系下，当底部边界条件固定后，经过充分的时间松弛可以达到稳态分布.值得注意的是本文并没有像其他模型那样假设内边界上总压或者热压为常数来约束温度或者密度，而是利用了二者与速度之间的相关关系.尽管如此，总压在内边界上的梯度比较小，在靠近边界的地方不足以驱动等离子出现大规模非径向流动，这与非径向速度和磁场的设定是相容的，而且模拟结果也表明，这种方法更能体现密度和温度的真实分布.

3 模拟结果

利用上述模型，首先模拟了2007年(CR2051-CR2065)的背景太阳风.这一时期处于第23太阳活动周的低年，近地观测到的日冕物质抛射事件较少(Richardson and Cane, 2010；Chi et al., 2016)，太阳风的结构主要受高速流调制，因此适合模型的测试.边界条件上的可调参数的具体取值在表 2中列出，除了WSA关系中的a₂在CR2051-CR2060与CR2061-CR2065采用不同的取值, 其余参数取值均相同.作为示例，图 1展示了CR2053的径向速度V_r、数密度N、温度T以及径向磁场B_r在内边界上的分布.在这一时期，电流片主要分布在南北纬30°以内，黄道面内为典型的四扇区结构，高纬度区域主要分布着源自极区冕洞的高速流，而在低纬度地区高速流和低速流相间分布.

图 2展示了CR2053太阳风传播到1AU之后在卡林顿坐标系中的分布，可以看到太阳风的物理量与内边界相比发生明显变化，最直观的是各物理量在经度方向上有向左大约50°偏移，这是共转效应造成的.此外还有一些其他变化：(1)太阳风速度平均有50 km·s^-1左右的增幅，这是因为在行星际随着日心距离的增加，磁场和温度逐渐衰减，大部分的磁能和内能转换成了太阳风的动能，而WSA经验模式无法体现这一过程，正因如此模型在内边界上的速度在WSA关系的基础上再减去50 km·s^-1；(2)在低纬度区域，由于高速流与低速流之间的相互作用形成数个压缩区和稀疏区(Gosling et al., 1972; Gosling and Pizzo, 1999)，在压缩区密度的局部值超过30 cm^-3，磁场强度达到15 nT以上，远高于各自的平均值；(3)尽管温度和速度仍然保持正相关关系，但是在低纬度区域高速流前缘压缩区的温度大于高速流后缘稀疏区的温度，表明在压缩区存在显著的压缩加热效应.

为了检验模型的可靠性，将模拟值和近地卫星观测进行了对比.太阳风观测数据从OMNI网站(http://omniweb.gsfc.nasa.gov/)下载而来，它包含了WIND和ACE飞船在L1点的实地观测.模拟中经度和纬度方向上的网格的分辨率为1°×1°，对应的时间分辨率~1.8小时.当卡林顿周开始时，L1点位于经度360°处，然后逐渐减小，最终在卡林顿周结束时到达0°位置.考虑黄道面和赤道面的夹角，通过插值得到L1点每小时的模拟值，然后和每小时平均的观测值对比.图 3展示了模拟结果和观测的时间序列的对比图，从上到下依次为径向速度V_r、数密度N、温度T以及磁场强度B，其中红色线代表模拟值，蓝色线代表观测值.由于模型中的网格分辨率有限，并且数值解只能识别若干个网格以上的空间尺度的结构，所以模拟的曲线显得比较平滑.这一时期捕捉到几乎所有的高速流，而且高速流的峰值和持续时间与观测比较一致.对于空间天气预报，高速流的到达时间也是非常重要的.根据交叉相关分析，模拟的全年速度的时间延迟为0小时，也就是说虽然部分高速流的到达时间与观测不同步，但是模拟结果不存在系统的时间延迟.注意到当速度和观测吻合较好时，尤其当高速流的到达时间与观测同步时，密度和磁场的压缩结构与观测符合也比较好，反则反之，这也说明了内边界上速度分布的重要作用.

表 3给出了这一时期模拟和观测的一些统计分析结果，、模拟平均值、均方根误差、平均相对误差以及相关系数.可以看到，所列的四个参数模拟和观测的平均值都非常接近.速度的均方根误差为107 km·s^-1，平均相对误差为18%，相关系数也达到0.60.数密度的均方根误差大致和观测值的平均值相当，相关系数也只有0.32，尽管如此，在高速流和低速流之间的变化幅度与观测基本一致.在内边界上，如果采用N∝V^-2，调节动量通量的取值使模拟结果与L1点观测符合较好时，极区的密度要比Ulysses观测低20%左右；如果采用N∝V^-1，调节质量通量的取值使模拟与L1点观测符合较好时，极区的密度要比Ulysses观测高75%左右.可以说，式(8)更好体现了0.1AU处密度和速度之间的关系.温度的均方根误差接近观测值的平均值，平均相对误差也达到了87%，为四个参数的最大值.磁场强度的模拟值与观测的平均值相差5%，这表明着内边界上磁场强度的设置是比较合理.同一时期用Gong的磁图作为PFSS+SCS模型的输入反演得到的磁场仅为80 nT左右，也就是低估了大约75%的磁场通量.以CR2053模型所确定的边界条件为例，磁压为动压的30%到40%，热压为动压的7.5%，可以说磁场在行星际动力学过程中(至少在内边界附近)地位是不容忽视的，磁压的作用一般是要强于热压的.如果磁场的强度低估两倍以上，那么模型所描述的状态将是与真实情况有明显差别.

Gressl等(2014)分别利用MAS/MAS、MAS/ENLIL以及WSA/ENLIL三种模型模拟了同一时期的背景太阳风.其中MAS/MAS和MAS/ENLIL明显高估低速流的密度，而WSA+ENLIL则低估了温度.三种模拟所得的磁场强度都至少与观测相差两倍.与这些模型相比，本文的模型通过调节非常少的参数，所得的密度、温度以及磁场强度与ACE/WIND观测的平均值、变化幅度都比较一致，而且各个参数的相关系数要高0.1左右，因而更能反映背景太阳风在行星际空间的三维分布.

为了进一步检验模型，本文继续模拟了2016年(CR2171-CR2185)的背景太阳风.这一时期位于第24太阳活动周的下降期，行星际磁场和密度都相比2007年有了显著的增加.作为示例，图 4和图 5分别展示了CR2180太阳风的四个参数在0.1AU和1AU的分布.从图中可以看到这时期电流片延伸到高纬度地区，黄道面的位置只有两个扇区结构；速度的分布比较复杂，无论是在高纬度还是在低纬度区域，都有高速流和低速流分布，而且这些高速流源自许多大小不同的冕洞.在1AU, 共转相互作用导致的稀疏区和压缩区的分布也扩展到中纬度区域.

图 6展示了2016年的1AU处模拟值和观测值的对比.前半年有7个高速流未被模型捕捉到，相应的密度和磁场的压缩结构也没有出现，这再次说明了速度分布的重要意义；后半年四个参数的变化幅度和峰值也与观测很好相符.这一年的分析结果如表 4所示，无论是相关系数还是均方根误差都与2007年大体一致.可以说，通过调节个别参数，该模型可以重现不同太阳活动周时期的太阳风的空间结构.

4 总结

本文建立了一套三维MHD数值模型，该模型的内边界设置在0.1AU处，在六片网格系统下利用TVD Lax-Friedrich格式求解理想MHD方程组，从而模拟背景太阳风在行星际空间的传播和演化.模型以GONG的观测磁图作为输入数据，根据PFSS+SCS模型得到内边界上的磁场分布，采用WSA经验模型得到内边界上的速度分布，结合卫星观测来确定密度、温度的分布，从而构成了一组完备自洽的边界条件.在边界条件中保留了6个可调参数，以便模拟不同太阳活动期的太阳风，并尽可能使各个参数的变化特征贴近观测.经过测试，模型可以高效稳定运行.本文利用该模型分别模拟了2007年和2016年的背景太阳风，模拟结果在1AU不仅捕捉到了大多数的高速流，而且密度、温度和磁场强度也没有明显高估或者低估.相比与WSA+ENLIL等模型，该模型能更好反映背景太阳风在行星际空间的三维分布.

值得注意的是，目前使用单幅磁图代表整个卡林顿周期内的状态，而实际上磁图是在不断变化的，在个别卡林顿周期的起始时刻，例如CR2058(2007年6月21日)的开始，各物理量的时间序列图上有明显的间断.为了更加准确预报太阳风参数的变化，未来将考虑采用每小时更新的GONG的磁图驱动模型，建立时变的三维数值模型，模拟更加接近实际的背景太阳风，并在此基础上研究扰动事件的传播过程.

致谢  我们使用了从GSFC/SPDF OMNI Web界面http://omniweb.gsfc.nasa.gov获取的太阳风数据，从NSO GONG网站http://gong.nso.edu/获取的太阳磁场概略图，以及Richardson和Cane的近地ICMEs列表http://www.srl.caltech.edu/ACE/ASC/DATA/level3/icmetable2.htm，对此表示感谢.