行星热演化，即行星内部温度和热流的随时间变化，从能量源头上控制了行星内部物质运动和动力学演化，因而是行星物理学研究的重要内容。在类地行星内部，自从早期核幔分异完成之后，核幔边界就成为最主要的化学成分界面，行星核与行星幔之间不再进行物质交换，但核幔之间的能量交换则一直持续不断。由于缺乏对行星核幔边界温度和热流的有效观测约束，行星热演化无法分成行星核与行星幔2个独立系统分别加以研究。换言之，行星热演化问题不能核幔解耦，必须进行核幔耦合的研究。

迄今为止，文献中发表的行星热演化的核幔耦合研究成果还为数不多。究其原因，主要是受制于缺乏适合的数值模型。一方面，基于精确模型的直接数值模拟计算量过大。例如，地球、金星等行星核与行星幔内部均存在强烈的对流^[1-2]，行星热演化的核幔耦合直接数值模拟对计算资源要求过高，目前还无法实现。另一方面，以过于简化模型为基础的研究难以充分反映内在物理而具有严重的局限性。已有的行星热演化的核幔耦合研究往往采用集中质量点模型^[3-5]，仅考虑行星核与行星幔的平均温度演化，忽略了核幔内部温度的空间变化及其与热流间的内在联系。本研究以水星作为研究对象，探索繁简适度、计算可行的行星热演化数值模型。

水星是太阳系八大行星中最小的类地行星，大约形成于40亿年前，密度与地球相近，表面大气稀疏，内部结构以幔薄、核大(核半径约为行星半径的80 %)为鲜明特征^[6]。水星独特的表面条件和核幔结构特点使之成为研究行星热演化的理想对象。一方面，由于水星大气稀疏，其表面温度可由行星辐射平衡温度较好地近似，因而表面边界条件简单。另一方面，水星幔层很薄，其内部热对流很可能处于层流而非湍流状态，易于合理近似。

近年来，水星探测器水手10号(Mariner 10)和信使号(MESSENGER)，均探测到水星表面存在磁场。然而，观测到的水星磁场异常微弱，其强度仅为地球现今磁场的百分之一，是行星发电机理论在科里奥利力与洛伦兹力平衡假设下预测的水星磁场的1/30^[7]。对水星表面磁场异常微弱的一种解释是，核幔边界附近的水星核区域处于亚绝热状态，不发生对流，因而对水星磁场没有贡献^[7-8]。这样的解释虽非学术界定论，但从中可以清楚地看出行星核的结构与所处热状态对行星磁场强度的控制作用。对水星表面磁场异常微弱的另一种解释是，与自内向外进行的地核结晶过程不同，水星核凝固是自外向内进行的。近年来基于不同行星的物性参数比较估计，一些学者提出了水星核自外向内凝固的可能性^[9-10]。水星核不同的凝固方式会对水星磁场具有决定性的影响。Stanley等^[11]通过假设水星核内凝固层的存在，得到的水星发电机数值模拟结果初步解释了水星磁场特征。此后，众多水星磁场模型都表明，水星核中凝固层产生的磁屏蔽效应对磁场特征影响至关重要^[7]。因此，水星核凝固过程是水星磁场研究的基础，能够为水星发电机模型提供重要的约束。本研究基于核幔耦合的数值模型，全面比较和分析水星核向外凝固与向内凝固两种模式下的水星热演化过程，为揭示水星磁场的成因提供启发和参考。

1 水星核凝固模式

一般地，随着类地行星核的逐渐冷却，在温度降到熔点(凝固点)的地方将首先发生从液态到固态的结晶凝固。行星核凝固方向是从内向外还是从外向内，取决于是行星中心还是核幔边界先冷却到相应位置的熔点温度。在行星核中，受压力控制的熔点总是从行星中心向外逐渐降低，而固相线梯度可定量地由Lindeman公式^[12]给出

$ \frac{\partial T_{\mathrm{L}}}{\partial r}=-2 \lambda r T_{\mathrm{L}}, $

(1)

$ \lambda=\frac{4 \mathsf{π} G \rho_{\mathrm{c}}^2(\gamma-1 / 3)}{3 K}, $

(2)

其中：r为距行星中心距离，T_L、ρ_c、γ、K分别为行星核的熔点温度、密度、Grüneisen参数和体积模量，而G为万有引力常数。

类地行星的液核由于黏度非常低，总是处于激烈的对流状态，其内部温度的径向变化可由绝热线近似(忽略非常薄的液核边界层)，对应于绝热温度梯度^[12]

$ \frac{\partial T_1}{\partial r}=-2 \nu r T_1, $

(3)

$ \nu=\frac{2 \mathsf{π} G \alpha_{\mathrm{c}} \rho_{\mathrm{c}}}{3 c_{p, \mathrm{c}}}, $

(4)

其中: T_l为液核温度，α_c与c_{p, c}分别为液核热膨胀系数与定压比热。

由于液核绝热梯度与固相线梯度具有相似的对压力依赖关系，而且在相变发生时T_l=T_L，行星核凝固方向是从内向外还是从外向内，可通过绝热梯度与固相线梯度的比值

$ \frac{\partial T_1 / \partial r}{\partial T_{\mathrm{L}} / \partial r}=\frac{\nu}{\lambda} $

(5)

小于1或大于1来判断^[13]。具体而言，若ν < λ，则固相线是超绝热的，具有比绝热线更陡的径向变化趋势，因而随着液核逐渐冷却，近似沿绝热线变化的液核温度剖面一定首先在行星中心与固相线相交，导致行星核凝固自内向外发生^[13]。相应地，行星核在凝固过程中由固态内核与液态外核构成，与现今地球类似。反之，若ν>λ，则固相线是亚绝热的，凝固将自外向内进行^[13]，凝固过程中的行星核由固态外核与液态内核构成。图 1示意性地绘出了本文研究的2种水星凝固模式。图中，a、b、c分别表示水星的行星半径、核幔边界(CMB)半径与内核边界(ICB)半径，T_s、T_l、T_m分别表示水星固态核、液态核与水星幔温度，而T_L为ICB所处的固相线温度。

图 1所示2种水星核凝固模式，在前人研究中均有涉及。Edgington等^[10]基于分子动力学模拟的结果，认为Fe₈₀S₁₀Si₁₀(下标表示物质的量所占百分比，即摩尔百分比) 组分下的水星液核顶部可能存在外凝固层，而Hauck等^[14]基于蒙特卡罗方法的研究结果认为，水星液核顶部可能存在固态富FeS层。这些研究都对应于水星核自外向内凝固的模式。Glassmeier等^[15]则考虑了水星核自内向外的凝固模式，认为水星磁场特征主要是基于水星磁层感生磁场与水星发电机磁场的相互作用产生的。此外，Tian等^[16]提出一种水星液核内外两侧均存在凝固层的发电机模型，并讨论了外凝固层以电磁屏蔽减弱水星磁场的作用。

需要说明的是，在前人研究中，仅从静态角度考虑水星核凝固的不同模式，重点在于研究水星核凝固模式与水星磁场特征间的关系。与前人研究不同，本研究系统分析在水星核幔耦合热演化的动态过程之中，水星核凝固的不同模式(图 1)如何控制水星内部温度与热流的演化，从而为研究水星磁场的产生与演化提供参考性基础数据。

2 水星热演化模型 2.1 水星幔传热计算

由于水星表面几乎没有大气，本文将水星表面温度设为恒定的辐射平衡温度

$ \left.T_{\mathrm{m}}\right|_{r=a}=T_0=445 \mathrm{~K}, $

(6)

其中: T_m为水星幔温度，a为水星半径。

由于水星核冷却和凝固所释放的热流持续对水星幔底部进行加热，而放射性元素衰变不断地在水星幔内部生热，水星幔中存在静止盖层对流^{[9-10, 18]}形态的热对流。由于水星幔黏度对温度的高度敏感性，低温高黏的盖层部分变形微弱，而对流运动主要发生在静止盖层之下的高温低黏区域^[9-10]。静止盖层对流的激烈程度可由Rayleigh数(Ra)

$ R a_{\mathrm{i}}=\frac{\rho_{\mathrm{m}} g_0 \alpha_{\mathrm{m}}\left(T_{\mathrm{CMB}}-T_0\right)(a-b)^3}{\kappa_{\mathrm{m}} \eta_{\mathrm{i}}} $

(7)

定量表征^[9-10]。其中，ρ_m、α_m、μ_m、κ_m分别为水星幔密度、热膨胀系数、黏度与热扩散系数，a、b分别表示水星的行星半径与核幔边界半径，T₀与T_CMB分别为水星表面与核幔边界温度，g₀为重力加速度。Ra定义式(7)中的η_i为水星幔底部黏度^[19]，其取值依赖于核幔边界温度

$ \eta_{\mathrm{i}}=\eta_0 \exp \left[-E\left(T_{\text {ref }}-T_{\mathrm{CMB}}\right)\right], $

(8)

式中, η₀为参考温度T_ref下的水星幔黏度，而黏度对温度的依赖性由参数E定量刻画。

热对流相对于热传导的强化传热效应可由Nusselt数(Nu)

$ N u_{\mathrm{m}}=\frac{Q}{Q_{\mathrm{c}}} $

(9)

定量表征。其中，Q和Q_c分别代表施加同样温差(T_CMB-T₀)下的对流热流与传导热流。为参数化表示水星幔对流传热Nu对Ra和跨水星幔温差的依赖性，本研究采用如下的定量关系^[20]

$ N u_{\mathrm{m}}=1.89 \theta^{-1.02} R a_{\mathrm{i}}^{0.2}, $

(10)

其中，θ为Frank-Kamenetskii参数(F-K参数)^[20]

$ \theta=E\left(T_{\mathrm{CMB}}-T_0\right). $

(11)

虽然式(10)给出的定量关系最初是从拟合底部加热矩形区域内二维静止盖层对流计算结果得到的，但因其与渐近边界层理论得出的结果高度一致^[20]，所以能够适用于更普遍的情形，包括内部加热情形与球壳几何。例如，内部加热内外径比0.55球壳内三维静止盖层对流的数值模拟结果与式(10)差别很小^[21]。所以，定量关系式(10)应用于内外径比高达0.8的水星幔应该是合理可行的。

式(9)与式(10)组合起来提供了在热演化模型中简化水星幔传热计算的一种可行方法。式(9)表明对流热流是同样温差下传导热流的Nu_m倍，而式(10)表明该Nu仅依赖于Ra和F-K参数。要计算水星幔热流，可以先求出给定温差和物性参数的Ra和F-K参数，再利用式(10)计算Nu，然后求解热导率放大Nu_m倍的球壳热传导方程

$ \frac{\partial T_{\mathrm{m}}}{\partial t}=N u_{\mathrm{m}} \kappa_{\mathrm{m}}\left(\frac{\partial^2 T_{\mathrm{m}}}{\partial r^2}+\frac{2}{r} \frac{\partial T_{\mathrm{m}}}{\partial r}\right)+\frac{H}{c_{p, \mathrm{~m}}}, $

(12)

$ \kappa_{\mathrm{m}}=\frac{k_{\mathrm{m}}}{\rho_{\mathrm{m}} c_{p, \mathrm{~m}}}, $

(13)

其中：k_m、c_{p, m}、H分别为水星幔热导率、定压比热与内部生热率。式(12)中，将水星幔模拟成热导率为Nu_m·k_m的材料，其对应的传导热流必然也放大Nu_m倍，恰好等于式(9)中对应的对流热流Q。这样，通过求解热传导方程(12)得到的热流与对流热流总体等效，从而不必直接数值模拟复杂的黏度依赖于温度的对流问题就可以分析水星热演化。

本研究考虑水星幔内随时间衰变的放射性生热，取水星核开始凝固时作为水星热演化的时间起点t=0，将内部生热率指定为

$ H=H_1 \mathrm{e}^{-\sigma\left(t-t_s\right)}, $

(14)

其中，H₁为水星幔内现今放射性生热率，σ为衰变常数，而t_s为水星固体核的现今年龄。

本研究采用的水星幔对流传热简化计算方法，实质上是将对流数值模拟的结果参数化为Nu依赖于Ra和F-K参数的定量关系，再与热演化模型结合。这为将文献中行星内部对流数值模拟研究的大量成果应用于构建行星热演化模型，开辟了广阔的空间。需要说明的是，对流强化传热是通过在流体中形成边界层结构实现的，利用式(9)~式(14)的计算，虽然忽略了对流的内部结构，但达到了与对流系统热流总体等效的效果，从而相对于直接数值模拟极大地简化了行星热演化的数值计算。

2.2 星核向外凝固的水星热演化模型

星核自内向外凝固的水星热演化模型，包括水星幔、液态外核和固态内核3个区域(图 1(a))。在水星幔内部，温度满足球壳热传导方程(12)，而在水星幔的外表面(即水星表面)指定恒温边界条件(6)。

在水星的内核边界r=c上，内外核温度连续，且等于熔点温度

$ \left.T_{\mathrm{s}}\right|_{r=c}=\left.T_1\right|_{r=c}=\left.T_{\mathrm{L}}\right|_{r=0} \cdot \mathrm{e}^{-\lambda c^2}, $

(15)

其中，$\left.T_{\mathrm{L}}\right|_{r=0}$为水星中心的熔点温度。在水星液态外核内部c < r < b，温度分布由绝热线近似为

$ T_1=\left.T_{\mathrm{L}}\right|_{r=0} \cdot \mathrm{e}^{v\left(c^2-r^2\right)-\lambda c^2} $

(16)

在水星核幔边界r=b上，核幔温度连续，从而有

$ \left.T_{\mathrm{m}}\right|_{r=b}=\left.T_1\right|_{r=b}=\left.T_{\mathrm{L}}\right|_{r=0} \cdot \mathrm{e}^{v\left(c^2-b^2\right)-\lambda c^2} . $

(17)

液态外核的能量守恒可表示为

$ \left.k_{\mathrm{s}} \frac{\partial T_{\mathrm{s}}}{\partial r}\right|_{r=c}-\left.N u_{\mathrm{m}} \frac{k_{\mathrm{m}} b^2}{c^2} \frac{\partial T_{\mathrm{m}}}{\partial r}\right|_{r=b}=q_{\mathrm{L}}+q_{\mathrm{G}}+q_{\mathrm{A}}, $

(18)

$ q_{\mathrm{L}}=\rho_{\mathrm{s}} L \frac{\mathrm{~d} c}{\mathrm{~d} t}, $

(19)

$ q_{\mathrm{G}}=\frac{2}{15} \mathsf{π} G \rho_{\mathrm{s}} \Delta \rho\left(\frac{3 b^5-5 b^3 c^2+2 c^5}{b^3-c^3}\right) \frac{\mathrm{d} c}{\mathrm{~d} t}, $

(20)

$ q_{\mathrm{A}}=-\frac{\rho_{\mathrm{s}} c_{p, \mathrm{c}}}{c^2} \int_c^b r^2 \frac{\partial T_1}{\partial t} \mathrm{~d} r $

(21)

这里q_L、q_G与q_A分别代表水星核凝固释放的相变潜热、引力势能与外核冷却释放的热能，而k_s、L、c_{p, c}、ρ_s分别为水星核热导率、相变潜热、定压比热与固态内核密度。与前人研究^[20]类似，本文只在引力势能q_G的计算中考虑液核凝固伴随的密度变化Δρ，而在其余计算中，内外核设定为同样的密度，即

$ \rho_{\mathrm{s}}=\rho_{\mathrm{c}} . $

(22)

水星核凝固开始于t=0时刻。此时，水星中心处于熔点温度$\left.T_{\mathrm{L}}\right|_{r=0}$，而液核内温度沿绝热线下降到核幔边界温度$\left.T_{\mathrm{L}}\right|_{r=0} \cdot \mathrm{e}^{-u b^2}$，再经过水星幔下降到表面温度445 K。为简化起见，本文假设初始时刻水星幔的冷却可以忽略不计，因而幔内初始温度T_m0的分布可由求解稳态热传导方程

$ \frac{\partial^2 T_{\mathrm{m} 0}}{\partial r^2}+\frac{2}{r} \frac{\partial T_{\mathrm{m} 0}}{\partial r}+\frac{H}{c_{p, \mathrm{~m}}}=0, $

(23)

唯一确定。进而可以计算在核幔边界处的水星幔初始温度梯度，再乘以放大Nu_m倍的水星幔热导率，就得到初始核幔边界热流。此热流与外核能量守恒公式(18)相结合，决定了内核的初始生长速率。

在如上给定的初始条件下，水星固态内核在t>0以后的生长过程定义了一个球体Stefan问题。在固态内核内部0 < r < c，温度满足球体热传导方程

$ \frac{\partial T_{\mathrm{s}}}{\partial t}=\kappa_{\mathrm{s}}\left(\frac{\partial^2 T_{\mathrm{s}}}{\partial r^2}+\frac{2}{r} \frac{\partial T_{\mathrm{s}}}{\partial r}\right) . $

(24)

在水星中心温度取极大值

$ \left.\frac{\partial T_{\mathrm{s}}}{\partial r}\right|_{r=0}=0. $

(25)

而内核边界则是一个随时间向外移动的边界，其温度逐渐降低，但一直维持在固相线上，直到液核完全凝固或内核内部发生热对流。本文采用有限差分法，数值求解与外核能量守恒以及水星幔传热完全耦合的球体Stefan问题，算法可靠性验证在2.4节给出。

2.3 星核向内凝固的水星热演化模型

星核自外向内凝固的水星热演化模型，包括水星幔、固态外核和液态内核3个区域(图 1(b))。在水星幔内部，温度满足球壳热传导方程(12)，而在水星幔的外表面(即水星表面)指定恒温边界条件(6)。

在水星固态外核内部c < r < b，温度满足球壳热传导方程

$ \frac{\partial T_{\mathrm{s}}}{\partial t}=N u_{\mathrm{s}} \kappa_{\mathrm{s}}\left(\frac{\partial^2 T_{\mathrm{s}}}{\partial r^2}+\frac{2}{r} \frac{\partial T_{\mathrm{s}}}{\partial r}\right), $

(26)

其中Nu_s为固态外核的Nusselt数，即在水星固态外核内存在对流时，以将材料热导率k_s放大Nu_s倍的方式等效描述通过固态外核的对流热流。

在水星核幔边界r=b上，温度和热流均连续，即

$ \left.T_{\mathrm{s}}\right|_{r=b}=\left.T_{\mathrm{m}}\right|_{r=b}, $

(27)

$ \left.N u_{\mathrm{m}} k_{\mathrm{m}} \frac{\partial T_{\mathrm{m}}}{\partial r}\right|_{r=b}=\left.N u_{\mathrm{s}} k_{\mathrm{s}} \frac{\partial T_{\mathrm{s}}}{\partial r}\right|_{r=b} . $

(28)

在水星固态外核的内表面，即内核边界r=c上，内外核温度连续，且等于熔点温度

$ \left.T_{\mathrm{s}}\right|_{r=c}=\left.T_1\right|_{r=c}=\left.T_{\mathrm{L}}\right|_{r=0} \cdot \mathrm{e}^{-\lambda c^2} . $

(29)

在水星液态内核内部r < c，温度分布由绝热线近似为

$ T_1=\left.T_{\mathrm{L}}\right|_{r=0} \cdot \mathrm{e}^{v\left(c^2-r^2\right)-\lambda c^2} . $

(30)

而液态内核的能量守恒表示为

$ -\left.N u_{\mathrm{s}} k_{\mathrm{s}} \frac{\partial T_{\mathrm{s}}}{\partial r}\right|_{r=c}=q_{\mathrm{L}}+q_{\mathrm{G}}+q_{\mathrm{A}}, $

(31)

$ q_{\mathrm{L}}=-\rho_{\mathrm{s}} L \frac{\mathrm{~d} c}{\mathrm{~d} t}, $

(32)

$ q_{\mathrm{G}}=\frac{4}{15} \mathsf{π} G \rho_{\mathrm{s}} \Delta \rho c^2 \frac{\mathrm{~d} c}{\mathrm{~d} t}, $

(33)

$ q_{\mathrm{A}}=-\frac{\rho_{\mathrm{s}} c_{p, \mathrm{c}}}{c^2} \int_0^c r^2 \frac{\partial T_1}{\partial t} \mathrm{~d} r $

(34)

星核向内凝固问题的初始条件与星核向外凝固问题类似，但内核边界是随时间向内移动的，这样固态外核的生长过程定义了一个球壳Stefan问题，同样可采用有限差分法求解。

2.4 Stefan问题差分解的验证

本文采用有限差分法离散球体或球壳Stefan问题，再使用显式格式TDMA算法^[22]求解离散后的代数方程。在本研究求解的Stefan问题中，移动边界的熔点温度随内核边界对应的压力变化而变化，因此属于变熔点Stefan问题。由于从文献中找不到变熔点Stefan问题的可靠解析解或数值解进行对比，本节利用常熔点球壳Stefan问题的Davis and Hill^[23]解析摄动解验证本研究差分解的精度和可靠性。图 2对比了2种方法在Stefan数^[23]Ste=10条件下得到的相变边界随时间移动的计算结果。2种方法给出几乎相同的结果，表明本文采用的差分解的可靠性。基于有限差分方法，不但可以直接将常熔点算法推广到变熔点算法，而且可以方便地实现水星核幔间的传热耦合。

2.5 计算模型设定与参数取值

本文主要对比分析水星核向外和向内2种凝固模式下的水星核幔耦合热演化。为此，结合水星资料[3, 6, 12, 17, 19, 24-30]选取模型参数及其取值。表 1列出所有模型通用的水星参数值。这些参数值完全确定了水星核的固相线，熔点温度从水星中心的2 515 K逐渐下降到核幔边界处的2 000 K，对应于固相线衰减参数λ=6.34×10^-14m^-2。

水星核向外和向内2种凝固模式分别对应于绝热线衰减参数ν小于或大于固相线衰减参数λ。参考文献[3, 24]，选取如表 2所示水星核热膨胀系数与比热的2组取值，分别定义水星热演化模型系AS和BS，对应的ν与λ之比分别为0.49和1.98。因此，模型系AS由自内向外凝固模型组成，而模型系BS由自外向内凝固模型组成。

因凝固方向不同，模型系AS与模型系BS对应于非常不同的初始温度剖面。模型系AS的初始温度从水星中心熔点温度2 515 K，沿比固相线更缓的绝热线逐渐下降到CMB初始温度2 247 K，而模型系BS的初始温度从CMB熔点温度2 000 K开始，在水星核内沿比固相线更陡的绝热线逐渐上升到水星中心的3 151 K。

本文研究的模型系AS与BS各包含3个模型：AS1与BS1模型只考虑水星核对水星幔的底部加热作用，不考虑水星幔内的放射性生热以及水星核凝固伴随的引力势能变化；在AS1与BS1模型基础上，AS2与BS2模型考虑水星幔内的放射性生热，但不考虑水星核凝固伴随的引力势能变化；在AS2与BS2模型基础上，AS3与BS3模型进一步考虑水星核凝固伴随的引力势能变化。这样的细分，可以有效地将水星幔放射性生热与水星核凝固伴随的引力势能变化对水星热演化的贡献独立出来。

本文的计算模型对水星幔区域与固态核区域各采用500个格点进行差分离散，而时间差分步长取作Δt=0.118 5 Ma。通过计算结果对差分网格尺寸和时间步长的敏感性检验，可以验证本文计算模型选取的网格尺寸和时间步长已足够小。如图 3所示，将计算模型的差分网格尺寸和时间步长再减半，已不会改变水星核向内凝固模型BS3关于内核边界ICB半径演化的计算结果。特别地，在t=4Ga时，本文计算模型得到ICB半径1 271.237 km，而网格尺寸和时间步长同时减半的计算得到ICB半径1 271.885 km。因此，BS3模型预测水星核凝固4Ga形成的固体外核厚628.763 km，其中包含的离散误差仅为千分之一量级，完全可以忽略不计。

3 计算结果与讨论

图 4绘出了6个热演化模型给出的水星内核边界半径随时间演化的计算结果。水星核向外凝固的AS系模型给出水星固态内核随时间不断生长，但内核半径生长曲线则逐渐变缓。水星核向内凝固的BS系模型给出水星液态内核随时间不断缩小，且内核半径衰减曲线逐渐变陡。图 4结果可概括为：无论是向外还是向内凝固，均有内核越小、凝固前缘移动越快的趋势。

图 4表明，水星幔内放射性生热显著降低了水星核的凝固速度。如向外凝固的模型AS1预测水星核经过3.6 Ga就能完全凝固，而模型AS2则预测即使经过4 Ga水星核顶部仍有约100 km尚未凝固。水星核向内凝固的模型结果显示了同样的趋势。水星核经过4 Ga凝固后，模型BS1预测液核半径已小于800 km，远小于模型BS2预测的约1 400 km液核半径。水星幔的放射性生热显著增加了对流水星幔的传热负荷，必然导致水星幔传输水星核输入能量的能力大幅下降。

与水星幔内放射性生热相比，水星核凝固伴随的引力势能Q_G变化对水星核凝固速度的影响幅度较小。Q_G对向内凝固的影响大于对向外凝固的影响，且具有相反的效应：在向外凝固的模型AS3中减缓凝固(图 4(a))，而在向内凝固的模型BS3中加快凝固(图 4(b))。这是因为向外凝固使得密度大的固态核物质更接近水星中心，从而释放引力势能，而向内凝固使得密度大的固态核物质更远离水星中心，从而增加引力势能，导致需要从水星核传出的热流减小，减轻了水星幔的传热负荷。

图 5绘出了不同热演化模型对应的Ra和Nu的演化。模型系AS给出水星幔Ra从初始值5.7×10⁹快速单调下降(图 5(a))，而模型系BS则给出水星幔Ra从2.4×10⁸开始的演化趋势显著受控于幔内放射性生热(图 5(c))，包含放射性生热的模型BS2与BS3均给出Ra先增后减的趋势。这说明水星幔的内部加热有效地维持了幔内对流的强度。图 5(b)与5(d)表明，模型系AS和BS对应的Nu均一直高于3，且具有与Ra相似的演化形态。虽然水星幔的静止盖层对流具有比热传导高2倍以上的传热能力，但总体来说，水星幔传热能力还是远低于包含板块构造的地幔对流。

图 6绘出了不同模型给出的水星中心、内核边界与核幔边界温度演化。对于水星核向外凝固的模型系AS，所有温度演化曲线都是单调下降的，表明水星内部经历持续冷却。如图 6(a)所示，模型AS1结果表明，在水星核完全凝固后，水星中心与核幔边界以更快的速度降温。这是由于星核凝固后，不再释放相变潜热，水星冷却随即加快了速度。

对于水星核向内凝固的模型系BS，水星中心温度随时间单调下降，说明液态内核在持续冷却，而向内逐渐移动的内核边界逐渐升温，则是熔点随压力升高的必然结果。图 6给出的核幔边界温度演化显示，模型BS1给出水星幔逐渐冷却，而模型BS2与BS3则给出水星幔先缓慢加热再逐渐冷却。这样的差异显然可由水星幔的放射性生热解释。

图 7对6个热演化模型给出了在t= 0、2、4 Ga时水星内部温度剖面及其与水星核内固相线的对比。本文计算表明，模型系BS给出的传导固态外核温度剖面总是对应于亚绝热的温度梯度。图 7(b)、7(d)与7(f)示出的温度剖面清楚地显示了从液态内核中的绝热线向固态外核中亚绝热温度梯度的转换。相应地，模型系BS的计算过程中，固态外核内部不会发生对流失稳，而总是以传导方式传热，对应于公式(28)中固态外核的Nu_s恒为1。

图 8给出不同热演化模型给出的通过水星表面、核幔边界与内核边界的热流演化，其中从ICB-到ICB+的热流增加表示水星核凝固释放的相变潜热，从ICB+到ICB++的热流变化表示水星核凝固伴随的引力势能变化。Surface-Rad表示扣除放射性生热之后的表面热流，其高出CMB热流的部分对应于水星幔冷却释放的热流。

图 8显示，在AS1与BS1模型中，初始地表热流均等于核幔边界热流，反映了初始水星幔处于无内部生热的稳态传热状态，而其余模型由于放射性生热的存在，初始时刻的表面热流比核幔边界热流高4 TW。模型系AS与模型BS1给出的曲线具有相同的热流顺序，即从ICB-、ICB+、CMB到水星表面热流依次递增，表明随着水星核凝固进程的发展，水星幔、外核、内核都在逐渐冷却，流出每个区域的热流均高于流入热流。然而，模型BS2与BS3则反映出伴随液态内核逐渐冷却，固态外核在一段时间内持续吸热，即流入热流高于流出热流。这与图 6和图 7给出的温度演化结果是一致的。

以下对比分析模型AS3与BS3结果对水星磁场的意义。根据行星发电机理论，行星磁场是由液核中的对流形成的，所以只有通过液核的热流才能对产生和维持水星磁场有贡献。向外凝固的AS3模型表明，在液核凝固40Ga的时间范围内，离开液态外核的CMB热流虽然随时间逐渐减小，但一直维持在1.8 TW以上。这样的热流水平虽然显著低于地球核幔边界现今热流的早期估计值3~4 TW^[31-33]和近年来估计值5~15 TW^[34]，但相差不到10倍。从外核热流与磁场的能量联系来看，有理由推断模型AS3给出的热流水平不能解释仅为现今地磁场强度1 % 的观测水星磁场。

与AS3模型不同，水星核向内凝固的BS3模型结果表明，离开液态外核的ICB-热流在水星幔内部放射性生热影响下呈现先增加后减小的演化，而固态金属外核则逐渐增厚，其对水星磁场的屏蔽作用也随之逐渐增强。文献[16]考虑水星核顶部稳定传导层的磁屏蔽效应，并通过重现观测水星磁场特征约束稳定传导层厚度为300~900 km。本文BS3模型给出的处于亚绝热状态的固态外核，不论是从传热效果上，还是从磁屏蔽效果上，均等同于文献[16]考虑的稳定传导层，因此，文献[16]估计的稳定传导层厚度可以看作是对现今水星固态外核厚度的定量约束。这一约束与BS3模型结果相结合，可以进一步约束水星固态外核年龄与流过液态内核的热流。对应于300~900 km的水星固态外核厚度，图 3显示水星固态外核年龄在2.8 Ga以上，而估计外核厚度的中间值600 km对应于外核年龄3.9 Ga。这是本文将水星幔放射性生热率公式中的t_s设成4 Ga的原因。图 8(f)显示，在水星核凝固2.8~4 Ga之间，对应的ICB-热流介于0.8 ~0.4 TW。如此低水平的液核热流与固态外核的磁屏蔽效应相结合，足以解释水星外部观测磁场微弱的特征。

4 结论

核幔耦合的行星热演化研究需要繁简适度、计算可行的数值模型。本研究利用Nu依赖于Ra和F-K参数的定量关系，对水星幔静止盖层对流传热进行热流等效的参数化近似计算，并以绝热线近似液态水星核内的温度分布，从而建立水星热演化的核幔耦合数值模型。在此基础上，全面比较和分析了不同星核凝固模式下的水星热演化过程。基于计算结果，得到如下认识和结论：

1) 水星核无论是向外还是向内凝固，均呈现内核越小、凝固前缘移动越快的趋势。

2) 水星幔内放射性生热显著降低了水星核的凝固速度。水星核凝固伴随的引力势能变化对水星核向内凝固的促进作用强于对水星核向外凝固的减缓作用。

3) 星核向外凝固的水星热演化模型AS3揭示，在水星核凝固的整个过程中，水星核幔边界热流虽然小于现今地球核幔边界热流，但相差不足10倍。这样的热流水平，不能解释远远弱于现今地磁场强度的实测水星磁场成因。

4) 利用文献[16]对水星核顶部稳定传导层的厚度估计，本文热演化模型BS3预测，水星现今固态外核年龄大于2.8 Ga，而流过现今水星液态内核的热流介于0.8~0.4 TW之间。这样低水平的液核热流与固态外核的磁屏蔽效应相结合，足以解释水星观测磁场的低强度特征。