2018, Vol. 61
2. 中国科学院计算地球动力学重点实验室, 北京 100049;
3. 中国地震局地球物理研究所, 北京 100081
2. Key Laboratory of Computational Geodynamics, Chinese Academy of Sciences, Beijing 100049, China;
3. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
地球磁场由液态外核中的磁流体动力学(Magnetohydrodynamics, MHD)过程所产生和维持，这一机制被称为地磁场发电机(Roberts, 2007).发电机数值模拟意在通过数值计算，得到自持式发电机过程(Buffett, 2000).从1950年代到1990年代，发电机理论已经证实了外核发电机过程的可行性和基本条件.随着MHD数值模拟技术的成熟与大规模高性能计算机的飞速发展，发电机模拟在1995年迎来了突破.Glatzmaier和Roberts(1995b)的发电机模型GR95，不仅得到了三维自持发电机解，也展现了地磁场的一些基本形态特征和物理过程，例如偶极占优，地磁西漂和磁极倒转.此外，GR95模型还研究了内外核黏性耦合以及电磁耦合给内核差动旋转带来的影响(Glatzmaier and Roberts, 1995a)，并被随后的地震波研究所证实(Song and Richards, 1996).其后陆续有新的发电机模型被建立(Kageyama et al., 1995; Kuang and Bloxham, 1997; Christensen et al., 1998; Fearn, 1998; Dormy et al., 2000; Jones et al., 2000; Busse, 2000; Zhang and Schubert, 2000; Kono and Roberts, 2002; Glatzmaier et al., 2002; Takahashi et al., 2005)等等.大多数模型可以近似吻合地磁场的空间谱、长期变化以及地磁场倒转等现象.近年来，科学家已经能够逐渐运用发电机数值模拟结果解释地球及行星磁场的一些特殊的性质与机制(Jones, 2011).
地磁场发电机的数值模拟还有很多问题亟待解决，其中一个重要的障碍就是模型参数与真实参数的差异问题(Christensen and Aubert, 2006).对于大多数发电机模型而言，存在着4个主要的无量纲参数：表征热驱动大小的Rayleigh数Ra，表征黏滞效应与旋转效应之比的Ekman数E，表征黏性扩散与热扩散效应之比的Prandtl数Pr和表征黏性扩散与磁扩散效应之比的磁Prandtl数Pm(Christensen et al., 1999).受限于计算能力，发电机模型难以求解快速旋转以及低黏性参数(相对应较小的Ekman数)，目前所取参数(E~10-5)与估计参数(E~10-15)甚至到达近10个量级的差异(Christensen, 2008).这一鸿沟依靠计算技术的发展短期内仍难以解决.Davies等(2011)估计即便模型所用黏性与估计黏性仅相差6个量级(E~10-9)，也需要54000个可用的进程运行至少13000天来得到一个磁扩散时间的数值，这对于数值计算来说是个巨大的挑战.这些参数取值对于发电机模拟十分重要，影响着地磁场产生和维持的基本物理过程.因此，研究方程中参数对系统的影响，得到参数变化与磁场、流场和温度场等变化之间的渐进关系，是我们所关注的主要问题.
针对这一问题，前人已经积累了一些研究，其中有关于某一个参数的，如Ra(Olson et al., 1999; Kuang et al., 2008; Wang et al., 2013; Sreenivasan et al., 2014)，E(Sarson et al., 1998)，和Pm(Simitev and Busse, 2005; Calkins et al., 2012)；也有一些是系统性研究各个无量纲参数(Christensen et al., 1999; Christensen and Aubert, 2006; Jones, 2011).研究得到了许多比例关系及定量结论，例如发电机磁场在什么参数区间表现为偶极场占优，何时由多极场主导(Jones, 2011)；如何由给定参数估算流场和磁场的强弱(Starchenko and Jones, 2002)，以及磁场强弱是否与各种扩散系数相关(Christensen and Aubert, 2006)等等.但是4个主要的无量纲参数中，往往包含着共同的物理系数，因此其变化并不独立，这就导致在理解研究结论时，我们必须注意其适用性.此外，给出的一些比例关系往往涉及到多个额外定义的无量纲参数，不仅表现形式复杂，也给理解其物理机制增加了难度.集中研究某一项物理系数的影响，将会对以上问题提供必要的补充.
以黏性ν为例，目前单独研究其影响的数值结果还不多见.由于ν被包含在Ekman数E的定义之中，有时也以对E的研究等价于对ν的研究.Sarson等(1998)研究了E的变化对发电机的影响，发现当E减小时，必须通过增加Ra来维持磁场，并且此时磁场和流场都集中在内外核边界(ICB)附近.近来King等(2013)通过综合理论推导及数值模拟结果给出了一个简单明确的比例关系：u~C1/2E1/3，其中C为对流能量.虽然结果很简洁，但这个比例关系是基于无量纲参数Reynolds数(Re)和Ekman数(E)的，并且C、Re和E的定义式中均含有黏性项ν，这些参数关于ν是协变的.所以其实这些研究并没有直接给出黏性项的定量影响，也并未对黏性影响的机制提供直观的物理解释.
为研究黏性变化对发电机数值模型的影响，并且推测当模型参数接近估计参数时的输出结果，我们使用了MoSST模型进行数值模拟.MoSST模型发展自KB97模型(Kuang and Bloxham, 1997; 1999)，除了数值算法的一些改进之外，其在KB97模型的基础上也增加了一些额外的模块功能.该模型也应用于火星(王天媛等，2009；2013)和木卫磁场(焦立果等，2011)的起源研究.大多数的Geodynamo模型方程组无量纲化所采用的时间尺度值为黏性扩散时间，而MoSST模型采用的时间尺度值为磁扩散时间，所以无量纲参数中只有Ekman数中包含黏性项.因此用MoSST模型来进行黏性影响研究的优势在于：研究Ekman数的变化可实质性地等价于研究黏性的变化，且不涉及到其他参数的协同变化.本文的第二部分给出模型介绍，以及无量纲化的具体方案.第三部分通过固定其他参数值，改变E的取值范围(10-7~10-4)，得到了模型输出结果，并考察了黏性对外核流场、磁场和温度场的影响.第四部分主要从Geodynamo模型的物理机制上分析了数值结果；最后一部分是总结及讨论.
1 MoSST模型介绍球坐标下发电机模型的控制方程组由动量方程(Navier-Stokes方程)、磁感应方程、热方程以及速度场和磁场的无散约束所构成，其中速度场采用了Boussinesq近似：
|
(1) |
式中u为外核内流体的流动速度，Ω为地球自转角速度，P为修正后的压力(
以外核半径ro为长度尺度，磁扩散时间τ=ro2/η为时间尺度(η为磁扩散系数，η=1/μ0σ, μ0为真空磁导率，σ为外核内导电流体的电导率)，B=(2Ωμρη)1/2为磁场尺度，hTro为温度尺度(hT为在内核边界上的温度梯度：
|
(2) |
其中，Ω=Ωlz, Θ为温度扰动值，T0为CMB(核幔边界)温度.式中的无量纲参数(Ro为磁Rossby数，qк为磁Prandtl数)定义如下(本文第三部分表 1给出更加详细说明)：
|
表 1 MoSST模型各无量纲参数定义 Table 1 Definitions of dimensionless parameters in MoSST |
|
(3) |
可见在4个无量纲参数定义式(3)中，ν仅出现在E的定义中，这显著区别于其它模型(Kono and Roberts, 2002)的无量纲化方案.Geodynamo的数值模拟主要是为了理解磁场的起源和变化，关注的是在磁场发生显著变化的时间尺度内的MHD过程，因此选用磁扩散时间作为典型时间尺度是合宜的.
模型中所选用的初边界条件如下：速度场在ICB和CMB均为自由滑移边界条件ln·u=ln×(σν·ln)=0，其中ln为边界法向量，σν为黏性应力张量；磁场在ICB和CMB均为有限电导边界条件[B]=[ln×J]=[ln×E]=0 ([]表示通过边界的区别，E为无量纲电场)，在CMB外面还包含了一个电导率为其1/10的D″层；温度场在ICB和CMB均为固定热流边界条件∂θ/∂r=0.模型的初始状态来自我们之前的模拟结果.
MoSST模型在半径方向采用4阶紧致有限差分法求解，使用Chebyschev多项式展开作为径向网格配置节点.对于内核，外核和D″层，选取的节点数分别为35，39和19.在球面上，采用球谐展开和拟谱法求解，用FFT进行快速谱变换.球谐展开的截断阶数分别为：lmax=33, mmax=21.根据截断阶数，设置球向网格节点数分别为：θ向50个，ϕ向64个.模拟中所求解的网格数一共为277830个.
为避免高阶球谐截断对系统的影响，以及减少从初始状态到最终状态之间的过渡时间，MoSST模型中引入了超扩散(Kuang and Bloxham, 1999)，定义如下：
|
(4) |
式中l为球谐阶数.同样的，引入热超扩散和磁超扩散.该处理方法由Glatzmaier and Roberts(1995a)得到首个三维自洽发电机数值模型(GR95)时所引入.但以上形式相对于GR95，已经是一种弱形式.本文模拟中设定l0=4, ε=0.05.
2 数值模拟结果Geodynamo数值模拟的计算结果都受控于无量纲参数的影响，MoSST模型有四个相应的无量纲参数：Ro，E，Ra和qκ，具体说明如表 1所示.
表 1中的外核估值根据参数定义计算得到，所用到的扩散系数等物理参数的估值可参考Olsen(2007)和Christensen和Wicht(2007)的研究结果.从表 1可以看出，这些参数在数值模拟中的取值与实际值仍存在着较大的差距，目前的计算能力还远远无法达到参数的估计数量级.本研究已尽可能采用接近估计值的参数，数值计算在天河1号上完成.
2.1 诊断输出由于发电机方程组是一个各物理场高度耦合的非线性系统，当初始状态输入以后，在各无量纲参数的约束下，系统需要一段时间进行演化迭代以达到与所选参数相对应的最终状态.判断系统是否已经演化到终态，依赖诊断输出而实现.图 1给出了一个诊断输出图的例子，横坐标为演化时间，以磁自由扩散时间τf=r02/π2η为度量尺度(取r0=3.5×106m, η=2 m2·s-1, τf≈2×104年)；图 1的纵坐标分别表示速度场(VF)、磁场(MF)和温度场扰动(TP)、以及三个物理场所对应的典型空间尺度:速度场典型尺度(VFLS)、磁场典型尺度(MFLS)、温度扰动典型尺度(TPLS).
|
图 1 诊断输出示例 Fig. 1 Diagnostic outputs example |
诊断输出中的物理量为整个外核内的均方根结果，以外核流速为例：
|
(5) |
其典型尺度也是均方根结果，由以下公式计算得到：
|
(6) |
结合图 1，我们判断稳定输出的根据是：所有物理量及其典型空间尺度都已稳定达到一个τf以上.最终取的输出结果为一个稳定的τf以上时间范围内(不包含开始的过渡阶段)的平均结果，由图中直线所示.本研究所采用的是串行计算，计算时间依参数选择、初始状态、超扩散大小以及硬件性能而定.一般而言，得到图 1所示的稳定输出结果需要4~5天时间，但在极端情形下，得到一组稳定输出需要两周甚至更长时间.
2.2 各物理场随Ekman数的变化通过改变无量纲参数E的大小，得到了一系列稳定的输出结果，如表 2所示.
|
|
表 2 不同Ekman数下的输出结果 Table 2 Outputs under different Ekman number |
从表 2可以看出，在本研究所选参数范围内，除Ekman数取最小值(E=5×10-7)外，都产生了明显的发电机过程.在发电机过程持续范围内(1.25×10-6≤E≤2.5×10-4)，流速u随E增加而显著增加，其幅度变化超过一个量级；除了E=5×10-5和1.25×10-4两处例外，磁场B随E增加而呈现出单调减小的趋势，但其变化幅度不大，最大不超过30%；温度扰动Θ随着E增加而单调增加，只是变化幅度很小(~6%).这表明即使是很小的温度扰动，就可以带来流速的显著增加；而流速的显著增加，并不必然带来磁场的显著增强，相反，磁场却出现了减弱的趋势.
对于几个物理量的典型空间尺度而言，在发电机过程持续的参数范围内，VFLS(lu)随E增加而单调增加(在E=2.5×10-5处例外)，变化幅度为53%；MFLS(lB)未随E增加而单调变化，在E=2.5×10-5时出现最高值0.0689，在E=2.5×10-6时出现最低值0.0640，整体变化幅度不大(~8%)；TPLS(lT)随E增加而单调减小，与Θ变化反向，变化幅度不超过16%.前文提到长度以外核半径ro来无量纲化，模拟中取ro=3500 km；模型中同时包含了一个半径为ri=1200 km的内核，因此流场的典型尺度相当于在外核内有2~3个波长；磁场在外核中约有10个波长，尺度约为流场的1/3~1/5；而温度扰动在外核中只有约一个半波长，比对流的典型尺度更大.几种物理场典型尺度基本上都是以大尺度为主，温度扰动的空间尺度要大于对流尺度，而对流尺度又大于磁场尺度.这在一定程度上反映出温度扰动激发对流，而对流又激发磁场的动力学过程.关于不同尺度下的动力学过程，以及能量在各种尺度之间的层叠，可参考Huguet和Amit(2012)和Calkins等(2015)的研究结果.
为更加直观的对表 2进行分析，将数据绘制成图，如图 2所示.从图中可以更清晰的看出黏性的变化对于Geodynamo的影响主要体现在流场上：黏性越强，对流越活跃.为跟相关研究进行对比，对变化曲线进行了直线拟合，结果如图中黑色实线所示.从拟合结果可以看出，流速和Ekman数(黏性)之间满足指数关系式：
|
图 2 各物理场值及其典型尺度随Ekman数的变化 实线为物理量，虚线为相应典型尺度，不同形状代表不同的物理场，如图右上角所示；速度场的线性拟合结果用三角形直线表示. Fig. 2 Variations of the physical fields and their typical length scales with Ekman numbers Solid lines are physical fields; Dashed lines are the corresponding typical scales; Different type of lines represent different physical fields which are shown on the up right corner; the line with triangles represent result of the linear fitting of velocity fields. |
|
(7) |
关于VFLS(lu)与Ekman数之间的关系，我们也进行了拟合，得到(未在图 2中标出)：
|
(8) |
外核流体运动受各种作用力的影响，作用力之间的平衡机制决定着外核内的动力学过程，也关系到对模拟结果的解释.接下来我们对不同黏性下的各种作用力进行分析.将方程组(1)中的动量方程进行标注，得到
|
(9) |
式中包含五种作用力：惯性力FI，科氏力FC，洛伦兹力FL，黏性力FV和热浮力FA(阿基米德力).由于在求解过程中首先对动量方程求取旋度后再进行计算，压力梯度项消失，这里暂不予考虑.
通过上一节计算得到的各种平均物理场及其典型空间尺度，可以对各作用力大小进行估算，例如黏性力的大小可估算如下：
|
(10) |
类似地，空间求导都以典型空间尺度(根据求导对象来选择对应的尺度)的倒数来计，忽略点乘和叉乘.由于方程最后达到准稳定状态，忽略惯性项中的时变部分.根据以上规则得到的各作用力变化，在图 3中给出.图 2给出了u变化的拟合结果，这里同样，也对几个变化显著的作用力：惯性力，科氏力和黏滞力进行了线性拟合，其结果由图中虚线和公式标注.
|
图 3 各作用力随Ekman数的变化 不同形状的线条所对应的力如图右下角所示，线性拟合结果由公式和虚线表示. Fig. 3 Force balance under different Ekman number Lines with different shape represent different forces which are showed on the down right corner. Linearly fitted results are shown by formulas and dashed lines. |
由图 3可知，在Ekman数变化的整个范围内，由于没有改变输入能量Ra，热浮力只是随着温度扰动而出现了少许增加(6%).洛伦兹力随着磁场波动而发生波动变化，除了E=5×10-5和1.25×10-4两处例外，FL随E增加而呈现出单调减小的趋势.由于FL是B2与lB的叠加结果，其变化幅度超过B的变化，达到了50%.变化最明显的是另外三种力，我们一一分析.
科氏力在黏性增加的过程中单调增加，拟合结果显示其指数因子为0.49，这与u的变化完全相同，因为科氏力项的估算仅取决于u.惯性力随ν的变化曲线更陡，其指数因子为0.78.这是因为类似于FL，FI的变化是由于u2和lu的共同变化所影响的.变化最明显的是黏滞力项，随着ν的增加，FV迅速增长，拟合指数因子达到了1.66.究其原因，FV中不仅包含了ν本身，还受u2和lu的影响.
在E变化过程中，科氏力一直占据主导地位.在E较小时，热浮力主要与之平衡；但是随着科氏力和惯性力的快速增加，最后FI达到了与FA相当的幅度，二者共同平衡FC的影响.黏滞力虽然出现了快速增长，但最后仍与FC相差两个半量级以上，与FI和FA也相差一个半量级.洛伦兹力是惟一随着E增加而减小的力，且在整个过程中并未发挥显著影响.这表明系统始终处于弱场发电机过程.
2.4 流场和磁场形态为了进一步分析黏性变化对物理场形态的具体影响，我们绘制了不同Ekman数对应的流场形态图，如图 4所示.图中给出的是子午面剖视图，其中a, b, c分别对应的是不同Ekman数下的输出，依次为E=2.5×10-6，5×10-6，1.25×10-5.每幅剖视图包含左右两个部分，其中左半边的颜色代表轴对称(m=0)环型流场，右半边的线条代表轴对称极型流场.每幅图都是一段时间内的输出平均后的结果.环型流场(uT)和极型流场(uP)由下式定义：
|
图 4 轴对称流场形态随着Ekman数的变化 Fig. 4 Variations of axisymmetric velocity field with the Ekman number |
|
(11) |
剖视图的左半边是速度的φ向分量uφ，右半边是流线图(Ψu)，定义如下：
|
(12) |
式中s=r·sinθ，为场点到地球自转轴的距离.
由图 4可知，随着E的增加，流速出现了显著增加；当E=2.5×10-6时，流动主要(冷色部分)集中在正切圆柱(轴向与自转轴重合，与内核在ICB赤道面内切，与CMB在近极区外接的柱体)侧面和上下两个底面，而当E逐渐增加时，流动向正切圆柱内集中，同时在赤道面向CMB扩展.这两个特征都体现在环型流场部分(左半边).对于极型流场(右半边)，不论是在正切圆柱内，还是在正切圆柱外，都可以看出流动尺度的显著增加.
参照图 4，也给出了轴对称磁场的分布图，见图 5.图a, b, c分别对应的是不同Ekman数下的输出，依次为E=2.5×10-6，5×10-6，1.25×10-5.从环型磁场来看，场强随着E增加而减弱，并有逐渐向正切圆柱内集中的趋势.从极型磁场来看，三者并无显著差异，但都可以看出CMB以外的磁场(只有极型场存在)都是偶极占优.
|
图 5 轴对称磁场形态随着Ekman数的变化 Fig. 5 Variations of axisymmetric magnetic field with the Ekman number |
先前对边界条件影响的研究(Kuang and Bloxham, 1997; Kuang, 1999)表明，边界上的黏性耦合对geodynamo过程影响显著：当忽略黏性耦合时，该过程发生在整个外核；而当在边界上施加强黏性耦合时，磁场产生的区域则集中到正切圆柱以内.本研究选取的是无黏边条件，可在图 5中的环型磁场部分看出geodynamo过程发生在整个外核区域内.但是随着ν的增加，正切圆柱内成为geodynamo过程的主要场所.这表明在边界上施加黏性同在外核内增加黏性，对发电机过程产生了相似的影响.不同的一点是：当在外核内增加黏性时，流场在赤道面产生了一定的延展.发电机区域的缩小，或许正是导致磁场减弱的主要因素.
3 模拟结果的解释和讨论外核流速随着黏性增加而增加，这一点明显不同于一般的流体力学理论或者数值计算结果，但先前关于geodynamo类似的研究却支持同样的结论.本研究与先前研究有何异同，对geodynamo数值模拟在接近估计参数下的情形有何预测等，都将在这里进行讨论.
3.1 黏性对临界驱动力和流动尺度的影响地球外核流体运动强烈地受到地球自转效应的影响，旋转使得外核内的动力学过程发生了根本性的变化.主导地核运动过程的一阶项是地转平衡(这一点也可在图 3中科氏力始终占优看出来)，即在(9)式中只保留科氏力和压力梯度项：
|
(13) |
对上式两边取旋度，得到
|
(14) |
表明流动是二维的，沿旋转轴方向不变，即Taylor-Proudman定理(Taylor, 1963)，此时对流图案是平行于自转轴且绕自转轴转动的一系列平行圆柱形对流卷.但在封闭系统中，这种形式的流动既不能向外传递热量，也不能维持发电机过程.因此T-P约束必须被打破，使得控制方程的二阶项处于非地转状态.忽略惯性力和洛伦兹力，仅考虑黏性影响，对N-S方程求旋度，得到涡度方程，再对涡度方程取径向分量，得到
|
(15) |
式中ω是涡度，在式(6)中已经给出定义.由式(15)可知，在旋转不变的情况下，结合边界上ur=0，ν越大，u也就越大.实际上受速度其他分量、方程右侧涡度、以及其他作用力的影响，u并不与ν成简单正比关系.
黏性的增加打破了T-P约束，也表现在驱动系统对流的临界Rayleigh数RaC上.Chandrasekhar(1961)关于稳态旋转热对流的理论分析表明：
|
(16) |
由于Pr=ν/к，因此式(16)可推导为RaC~ν-1/3或者RaC~E-1/3；Roberts(1968)随后经过推导得到RaC=O(E-1/3)以及近来的RaC≈8.696E-1/3 (Roberts and King, 2013)；Kuang和Bloxham (1999)关于RaC与E的数值模拟结果表明：RaC≈10.1E-1/3.这些都说明，在旋转系统下随着黏性的增加，对流变得越来越容易发生.在相同的对流驱动力下，对流也越来越活跃.
一些研究也对对流卷的直径lconv进行了分析，发现lconv~E1/3H，其中H=ro-ri为球壳厚度(Chandrasekhar, 1961; Johnson and Constable, 1998; Roberts and King, 2013; Aurnou et al., 2015).图 6给出了对流卷直径lconv以及前文提到的正切圆柱的示意图.结合图 4，黏性增加不仅导致流速增加，也使得流动尺度变大.
|
图 6 正切圆柱及对流卷示意图 Fig. 6 Tangent cylinder and convection columns |
虽然关于E对RaC和lconv的研究结果比较多，但只有King和Buffett(2013)给出了关于E对u影响的定量结果：u∝C1/2E1/3.只是这一结果并不能简单理解为外核速度与黏性的关系，因为他们得到的结果是基于Re~C1/2E1/3，其中Re=ul/ν为Reynolds数，C=PH3/Aρ0ν3为对流能量.若假设P与A均与ν无关，经过推导，可得：u∝ν-3.这一推论显然与实际相悖.
因此本研究是目前为止惟一真正给出u同ν定量关系的工作，虽然这只是一个经验关系.据King和Buffett(2013)的研究，u同E之间的关系取决于N-S方程中的作用力平衡.他们的结论是基于VAC近似(FV, FA, FC三者主导)，若采用CIA近似，结果将变为u∝E1/5；而若采用MAC近似，结果为u∝E1/2.如果他们的结论适用于u与ν的关系，那么我们的模拟结果接近MAC近似.只是在前文作用力平衡的分析中，在E变化的整个区间，系统几乎都由热地转(AC)模式所主导；仅当E接近10-4时，系统才进入CIA模式.
这里关于lu随ν的变化(指数因子0.06)显然要弱于其他研究的结论(1/3)，这是因为对流卷是理想近似下的结果，实际上受各种作用力的影响对流会变得复杂.这一点可从图 7中赤道面上ur的分布看出来，此时对流卷在径向被拉长和扭转，已不再是理想的圆柱形，不能再用其直径来完全代表对流的尺度.当T-P约束被打破以后，对流已经变成准三维：虽然存在径向分量，其幅度也仅是uφ的2%(对比图 4).但是这2%的径向流动对于产生和维持geodynamo过程非常关键.离开了这2%的ur，热量就不能向外传递，流动也将退化成二维的.根据Cowling(1933)反发电机原理，二维流动将不能维持发电机过程.
|
图 7 赤道面上径向速度ur的分布(E=5×10-6) Fig. 7 Distribution of radial velocity on the equatorial plane (E=5×10-6) |
下面深入来看磁场的变化.King和Buffett(2013)的研究并未给出磁场随ν的变化关系，但前人已有很多相关研究.Starchenko和Jones(2002)认为，是驱动力的大小和作用力平衡而不是扩散系数决定着流速和磁场的大小；Christensen和Aubert(2006)假定扩散效应在控制方程中并不起主要作用，发现是可得能量(由流量Rayleigh数RaQ表征)而不是作用力平衡决定磁场大小；随后，Christensen(2010)对所有定义磁场大小的研究进行了总结，进一步指出磁场大小取决于可克服欧姆损耗的能量(F，见文中定义)，与旋转和电导无关，旋转只是控制是否产生偶极磁场；Jones(2011)也采用RaQ的概念给出了速度和磁场大小的估算.
综上，现在不再以传统的作用力平衡来估计流速和磁场的大小，而是转向能量输入的角度.能量输入中也不是全部能量，而是除了热传导外，特指驱动对流的这部分能量.从Christensen(2010)的研究看来，要得到更加准确的比例关系，驱动对流的这部分能量还要细化.扩散效应在决定流场和磁场的角色中似乎变得越来越不重要.
本研究中，在E超过两个量级的变化范围内(1.25×10-6≤E≤2.5×10-4)，磁场仅变化了30%，这部分支持了上述研究结论，即磁场随黏性变化不大.当E≤5×10-7时，磁场湮灭了，说明ν在某个临界点仍显著地影响着磁场的产生过程.此外在整个模拟中，流场随着ν的变化而发生了显著改变.本研究证明扩散效应的影响不可忽略，但在一定范围内可作适当假设.
3.3 超扩散的影响超扩散在2000年以后发展的一些模型中已经可以避免使用，这时需要适当提高模型的时空分辨率.为了提高计算效率，以及研究在低E下的输出，本研究中保留了超扩散.超扩散主要是人为抑制了小尺度对流的发展，对大尺度流动影响不显著，这也是造成各种物理场的典型空间尺度都是大尺度的原因之一.关于超扩散的影响，很多研究也进行了分析，例如Zhang和Jones(1997)发现超黏性显著增加了等效Ekman数，使对流的动力学过程的产生发生了变化；Grote等(2000)发现超扩散会影响磁场的产生机制和倒转过程.为更加明确地认识超扩散在MoSST模型中的影响，这里也进行了大致评估.
在模拟时，公式(4)中选取l0=4, ε=0.05，因此超扩散对于各物理场球谐展开中4阶以上的球谐项都有影响，并且随着空间尺度的减小(即球谐级数的增加)，影响逐渐增强.采用的最大球谐阶数为lmax=33，根据公式(4)可得其对应的超黏性为
|
(17) |
由此可见，超扩散对于计算结果的影响不容忽视，超扩散的引入使得高阶球谐项的黏度大幅度地增加.通过改变式(4)中的系数l0，进行了一些模拟，发现当超扩散减小时，带来了流速的降低和磁场的增加.这在一定程度上影响本研究结论，例如说在图 2中流场和磁场的拟合斜率实际上要更陡一些.在引用本研究结果时，需予以考虑.
3.4 对模型采用估计黏性时的预测按横坐标从右往左看图 3，可对模型采用接近估计黏性时的输出进行预测.当E越来越小时，系统由CIA平衡逐渐变化到CA平衡(热地转平衡)；当E进一步减小到逐渐接近估计值时，FI和FC继续减小而FL逐渐增加，最终系统将过渡到MAC平衡，此时系统将处于强场发电机状态，这也是前人所判断的geodynamo的实际运行状态(Taylor, 1963; Braginsky, 1967).但是FL随E的减小而增加的十分缓慢，这就导致模拟中真正实现MAC平衡较为困难，需要适当增加驱动力的大小.以本研究做线性外推，即便令E~10-15，FL仍将与FA保持有近一个量级的差异.导致这一点的原因来自所选参数空间的局限性(真实情形需增加Ra，同时减小qк)，以及模拟中采用了超扩散(使得图中E比真实情形要小).
需要说明的是，我们给出的力平衡只是一个平均的结果，其会随着流体位置的改变而变化.此外，这里给出的也只是一个大致的估计，因为无论是求导、点乘、叉乘，都是采用了近似或者忽略的策略.此外，热浮力项中，其大小也会随着场点距离地心的大小而显著不同.在做精确的作用力平衡计算时，必须考虑以上因素的影响.这里的分析仅给出了流体受力状况的一个基本的和整体的估计.
4 总结和讨论为研究黏性变化对地磁场发电机的影响，使用MoSST模型，通过与黏性对应的无量纲参数E近三个数量级，进行了数值模拟.使用MoSST模型的优势在于，在其无量纲化方案中，是以磁扩散而不是黏性扩散作为典型时间尺度，因而ν的变化仅包含在参数E中，E的变化可完全视为等效于ν的变化.在此之前的相关研究中，ν通常被包含在多个参数中，导致其结果的模糊性.
模拟结果显示，在所选参数范围内，外核流速对ν的变化是最敏感的，其拟合关系式为u~ν0.49.这打破了先前一些比例关系给出的结果.磁场和温度扰动都随E而大致呈现单调变化，但变化幅度不大：在E变化超过两个量级的过程中，二者变化幅度均不超过30%.对于几个物理量的典型空间尺度而言，速度场尺度的变化也是最大的，随E增加而单调增加了一半以上；磁场和温度场尺度的变化不大，都不超过16%.对作用力的分析表明，系统始终处于科氏力占优的状态，但随着E的减小，系统有从VAC过渡到CIA，最后到达MAC平衡的趋势.对流场和磁场形态的分析也进一步支持了以上结果.
流速随ν的增加而增加有违流体力学常识，在这里是受到了地球自转的影响.地球自转使得外核流动强烈受到T-P约束的影响，使得流动沿旋转轴不变；而黏性增加有助于打破这一约束，使得临界驱动力降低，最终带来了流速的增加.因改变扩散系数而带来的流速增加并未带来磁场的增强，因为场强主要是受到有效驱动能量的控制，与扩散系数的关系不大.然而，扩散系数在一定临界值下会显著影响磁场，导致磁场的消失.在进一步推导磁场和流场的比例规律时，需考虑本研究结果的影响.
本研究也有一些局限性需要提出，例如在E变化过程中，其他三个无量纲参数都未变化，它们的改变会定量影响这里得出的结论；作用力平衡分析只是整体和大致估计，需要精确计算；这里只是模拟了弱场发电机下的情形，对于强场发电机还需要进一步研究；在数值模拟过程中，为了提高计算效率而保留了超扩散，而超扩散又确实会对研究结果造成一定影响等等，这些方面将在下一步工作中完善.
致谢感谢天河一号提供计算资源.
Aurnou J M, Calkins M A, Cheng J S, et al.
2015. Rotating convective turbulence in Earth and planetary cores. Physics of the Earth and Planetary Interiors, 246: 52-71.
DOI:10.1016/j.pepi.2015.07.001 |
|
Braginsky S I.
1967. Magnetic waves in the Earth's core. Geomagnetism and Aeronomy, 7: 851-859.
|
|
Buffett B A.
2000. Earth's core and the geodynamo. Science, 288(5473): 2007-2012.
DOI:10.1126/science.288.5473.2007 |
|
Busse F H.
2000. Homogeneous dynamos in planetary cores and in the laboratory. Annual Review of Fluid Mechanics, 32(1): 383-408.
DOI:10.1146/annurev.fluid.32.1.383 |
|
Calkins M A, Aurnou J M, Eldredge J D, et al.
2012. The influence of fluid properties on the morphology of core turbulence and the geomagnetic field. Earth and Planetary Science Letters, 359-360: 55-60.
DOI:10.1016/j.epsl.2012.10.009 |
|
Calkins M A, Julien K, Tobias S M, et al.
2015. A multiscale dynamo model driven by quasi-geostrophic convection. Journal of Fluid Mechanics, 780: 143-166.
DOI:10.1017/jfm.2015.464 |
|
Chandrasekhar S. 1961.
Hydrodynamic and Hydromagnetic Stability. Oxford: Oxford University Press.
|
|
Christensen U, Olson P, Glatzmaier G A.
1998. A dynamo model interpretation of geomagnetic field structures. Geophysical Research Letters, 25(10): 1565-1568.
DOI:10.1029/98GL00911 |
|
Christensen U R, Olsen P, Glatzmaier G A.
1999. Numerical modelling of the geodynamo:A systematic parameter study. Geophysical Journal International, 138(2): 393-409.
DOI:10.1046/j.1365-246X.1999.00886.x |
|
Christensen U R, Aubert J.
2006. Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields. Geophysical Journal International, 166(1): 97-114.
DOI:10.1111/j.1365-246X.2006.03009.x |
|
Christensen U R, Wicht J. 2007. Numerical Dynamo Simulations. //Schubert G, ed. Treatise on Geophysics. Amsterdam: Elsevier, 245-282, doi: 10.1016/B978-044452748-6.00134-6.
|
|
Christensen U R.
2008. A sheet-metal geodynamo. Nature, 454(7208): 1058-1059.
DOI:10.1038/4541058a |
|
Christensen U R.
2010. Dynamo scaling laws and applications to the planets. Space Science Reviews, 152(1-4): 565-590.
DOI:10.1007/s11214-009-9553-2 |
|
Cowling T G.
1933. The magnetic field of sunspots. Monthly Notices of the Royal Astronomical Society, 94: 39-48.
DOI:10.1093/mnras/94.1.39 |
|
Davies C J, Gubbins D, Jimack P K.
2011. Scalability of pseudospectral methods for geodynamo simulations. Concurrency and Computayion Practice and Experience, 23(1): 38-56.
DOI:10.1002/cpe.1593 |
|
Dormy E, Valet J P, Courtillot V.
2000. Numerical models of the geodynamo and observational constraints. Geochemistry, Geophysics, Geosystems, 1(10).
|
|
Fearn D R.
1998. Hydromagnetic flow in planetary cores. Reports on Progress in Physics, 61(3): 175.
DOI:10.1088/0034-4885/61/3/001 |
|
Glatzmaier G A, Roberts P H.
1995a. A three-dimensional convective dynamo solution with rotating and finitely conducting inner core and mantle. Physics of the Earth and Planetary Interiors, 91(1-3): 63-75.
DOI:10.1016/0031-9201(95)03049-3 |
|
Glatzmaier G A, Roberts P H.
1995b. A three-dimensional self-consistent computer simulation of a geomagnetic field reversal. Nature, 377(6546): 203-209.
DOI:10.1038/377203a0 |
|
Glatzmaier G.A.
2002. Geodynamo simulations:How realistic are they?. Annual Review of Earth and Planetary Sciences, 30(1): 237-257.
DOI:10.1146/annurev.earth.30.091201.140817 |
|
Grote E, Busse F H, Tilgner A.
2000. Effects of hyperdiffusivities on dynamo simulations. Geophysical Research Letters, 27(13): 2001-2004.
DOI:10.1029/1999GL011155 |
|
Huguet L, Amit H.
2012. Magnetic energy transfer at the top of the Earth's core. Geophysical Journal International, 190(2): 856-870.
DOI:10.1111/j.1365-246X.2012.05542.x |
|
Jiao L G, Kuang W J, Ma S Z.
2011. Numerical simulations on origin of Galilean moons' magnetic anomalies. Science China Earth Sciences, 54(11): 1754-1760.
DOI:10.1007/s11430-011-4276-0 |
|
Johnson C L, Constable C G.
1998. Persistently anomalous Pacific geomagnetic fields. Geophysical Research Letters, 25(7): 1011-1014.
DOI:10.1029/98GL50666 |
|
Jones C A, Soward A M, Mussa A I, et al.
2000. The onset of thermal convection in a rapidly rotating sphere. Journal of Fluid Mechanics: 157-159.
|
|
Jones C A.
2011. Planetary magnetic fields and fluid dynamos. Annual Review of Fluid Mechanics, 43(1): 583-614.
DOI:10.1146/annurev-fluid-122109-160727 |
|
Kageyama A, Sato T.
1995. Computer simulation of a magnetohydrodynamic dynamo.Ⅱ. Physics of Plasmas, 2(5): 1421-1431.
DOI:10.1063/1.871485 |
|
King E M, Buffett B A.
2013. Flow speeds and length scales in geodynamo models:The role of viscosity. Earth and Planetary Science Letters, 371-372: 156-162.
DOI:10.1016/j.epsl.2013.04.001 |
|
Kono M, Roberts P H.
2002. Recent geodynamo simulations and observations of the geomagnetic field. Reviews of Geophysics, 40(4): 4-1.
|
|
Kuang W, Bloxham J.
1999. Numerical modeling of magnetohydrodynamic convection in a rapidly rotating spherical shell:weak and strong field dynamo action. Journal of Computational Physics, 153(1): 51-81.
DOI:10.1006/jcph.1999.6274 |
|
Kuang W, Jiang W, Wang T.
2008. Sudden termination of Martian dynamo?:Implications from subcritical dynamo simulations. Geophysical Research Letters, 35(14): L14204.
DOI:10.1029/2008GL034183 |
|
Kuang W J, Bloxham J.
1997. An Earh-like numerical dynamo model. Nature, 389(6649): 371-374.
DOI:10.1038/38712 |
|
Kuang W J.
1999. Force balances and convective state in the earth's core. Physics of the Earth and Planetary Interiors, 116(1-4): 65-79.
DOI:10.1016/S0031-9201(99)00116-8 |
|
Olsen P. 2007. Overview. //Schubert G, ed. Treatise on Geophysics. Amsterdam: Elsevier Science Publishers.
|
|
Olson P, Christensen U, Glatzmaier G A.
1999. Numerical modeling of the geodynamo:Mechanisms of field generation and equilibration. Journal of Geophysical Research:Solid Earth, 104(B5): 10383-10404.
DOI:10.1029/1999JB900013 |
|
Roberts P H.
1968. On the thermal instability of a rotating-fluid sphere containing heat sources. Philosophical Transactions of the Royal Society A:Mathematical, Physical and Engineering Sciences, 263(1136): 93-117.
DOI:10.1098/rsta.1968.0007 |
|
Roberts P H. 2007. Theory of the geodynamo. //Schubert G, ed. Treatise on Geophysics, Amsterdam: Elsevier, 67-105.
|
|
Roberts P H, King E M.
2013. On the genesis of the Earth's magnetism. Reports on Progress in Physics, 76(9): 096801.
DOI:10.1088/0034-4885/76/9/096801 |
|
Sarson G R, Jones C A, Longbottom A W.
1998. Convection driven geodynamo models of varying Ekman number. Geophysical & Astrophysical Fluid Dynamics, 88(1-4): 225-259.
DOI:10.1080/03091929808245475 |
|
Simitev R, Busse F H.
2005. Prandtl-number dependence of convection-driven dynamos in rotating spherical fluid shells. Journal of Fluid Mechanics, 532: 365-388.
DOI:10.1017/S0022112005004398 |
|
Song X D, Richards P.
1996. Seismological evidence for differential rotation of the Earth's inner core. Nature, 382(6588): 221-224.
DOI:10.1038/382221a0 |
|
Sreenivasan B, Sahoo S, Dhama G.
2014. The role of buoyancy in polarity reversals of the geodynamo. Geophysical Journal International, 199(3): 1698-1708.
DOI:10.1093/gji/ggu340 |
|
Starchenko S V, Jones C A.
2002. Typical velocities and magnetic field strengths in planetary interiors. Icarus, 157(2): 426-435.
DOI:10.1006/icar.2002.6842 |
|
Takahashi F., Matsushima M., Honkura Y.
2005. Simulations of a Quasi-Taylor State Geomagnetic Field Including Polarity Reversals on the Earth Simulator. Science, 309: 459-461.
DOI:10.1126/science.1111831 |
|
Taylor J B.
1963. The magneto-hydrodynamics of a rotating fluid and the Earth's dynamo problem. Proceedings of the Royal Society A:Mathematical, Physical and Engineering Sciences, 274(1357): 274-283.
DOI:10.1098/rspa.1963.0130 |
|
Wang T Y, Kuang W J, Ma S Z.
2009. Numerical simulation of Martian historical dynamo:Impact of the Rayleigh number on the dynamo state. Science in China Series D:Earth Sciences, 52(3): 402-410.
DOI:10.1007/s11430-009-0034-y |
|
Wang T Y, Wei Z G, Jiang W J, et al.
2013. Martian magnetic field properties before the termination of its core dynamo. Science China Earth Sciences, 56(8): 1452-1458.
DOI:10.1007/s11430-012-4510-4 |
|
Zhang K K, Jones C A.
1997. The effect of hyperviscosity on geodynamo models. Geophysical Research Letters, 24(22): 2869-2872.
DOI:10.1029/97GL02955 |
|
Zhang K K, Schubert G.
2000. Magnetohydrodynamics in rapidly rotating spherical systems. Annual review of fluid mechanics, 32(1): 409-443.
DOI:10.1146/annurev.fluid.32.1.409 |
|
焦立果, 匡伟佳, 马石庄.
2011. 伽利略卫星磁异常起源的数值模拟. 中国科学:地球科学, 41(12): 1822–1828.
|
|
王天媛, 匡伟佳, 马石庄.
2009. 数值模拟火星磁场古发电机:Rayleigh数对系统的影响. 中国科学:地球科学, 39(2): 157–165.
|
|
王天媛, WeiZ G, JiangW Y, 等.
2013. 数值模拟火星古发电机湮灭前磁场特征. 中国科学:地球科学, 43(1): 155–162.
|
|