欧罗巴是木星的第二个卫星, 是目前最有可能存在地外生物的星体 (Greenberg, 2008).为此, 欧罗巴成为空间探测计划的重点对象 (Prieto-Ballesteros et al., 2011).欧罗巴表面温度最低可达零下180 ℃(Allu Peddinti and McNamara, 2015), 同时其表面遭受着源于木星的高能粒子流的持续轰击, 生物很难在如此严酷的环境中生存 (Greenberg, 2005).人们不得不将希望寄予欧罗巴深部.欧罗巴表层物质为冰, 近年来人们认为冰层之下可能存在水 (Carr et al., 1998; Kivelson et al., 2000; Spohn and Schubert, 2003), 并且最新的研究认为欧罗巴水层的氧化还原过程与地球相似 (Vance et al., 2016).这些证据表明欧罗巴水环境有利于生物生存.另一方面, 如果冰层使得水与上表面隔绝, 这将大大降低生物存在的可能 (Gaidos et al., 1999).厚的冰层可能发生对流, 这有利于表面与水的连通; 而过薄的冰层不利于对流的发生, 这限制了表面与水的连通 (Gaidos et al., 1999).可见, 欧罗巴冰层厚度直接制约着生物存在的可能性.

欧罗巴冰层厚度的重要性类似于地球地壳厚度的重要性.欧罗巴冰层厚度控制着其温度结构、流变特性, 进而制约着冰体构造特征及演化过程 (Kargel et al., 2000).目前, 人们对于欧罗巴冰层的具体厚度尚未形成统一认识:有人认为小于1 km, 有人认为大于30 km (Billings and Kattenhorn, 2005).由于缺少直接观测, 人们一般通过间接方法来约束冰层厚度.间接方法的理论基础包括:板的挠曲, 陨石坑形态分析, 热-力学分析等 (Pappalardo et al., 1999; Billings and Kattenhorn, 2005).欧罗巴陨石坑分布范围广泛, 所以陨石坑约束之下得到的冰层厚度可能更具有普遍意义.

陨石撞击冰层形成陨石坑之后, 陨石坑深度将逐渐减小, 形态将逐渐平缓, 此过程为陨石坑的松弛过程.陨石坑深度减小到原始深度的1/e (e≈2.71828, 自然指数) 所需要的时间定义为松弛时间 (Thomas and Schubert, 1986).对于相同的陨石坑, 冰层厚度不同则松弛时间不同.反过来, 可以利用松弛过程约束冰层的厚度.Thomas和Shubert (1988)对木卫三陨石坑的松弛过程进行了有限元模拟.对于欧罗巴, 前人 (Buck et al., 2002; O′Brien et al., 2002) 基于松弛的方法估计的冰层厚度小于20 km.他们模型中没有考虑应力对冰黏度的影响 (牛顿流体)(Nimmo, 2004).然而, 实验室冰体蠕变实验研究表明, 冰的黏度是应力的函数 (非牛顿流体)(Goldsby and Kohlstedt, 2001).基于上述, 考虑应力对黏度的制约, Nimmo (2004)运用解析方法研究了不同尺度 (10 km, 100 km, 1000 km) 地形的松弛过程.研究结果认为欧罗巴冰层厚度小于10 km.我们注意到Nimmo的模型没有针对实际构造特征 (比如陨石坑) 进行分析, 更重要的是, 该文章在公式及数据引用方面存在疏漏.本文将改进这些不足.

基于有限单元法, 本文将对欧罗巴陨石坑的松弛过程进行动力学模拟.本文模型将:依据冰体实验数据, 将冰视为非牛顿流体; 修正前人文中存在的瑕疵; 针对不同厚度的冰体, 分别计算松弛时间.最终提出本文支持的欧罗巴冰体厚度.

本文有限元热-力耦合计算流程如图 1所示, 此计算流程与我们以前的冰体演化过程的有限元计算流程 (Yang and Shi, 2015) 相似.“初始化”中将建立计算几何模型、剖分网格、施加边界条件等.冰层的力学行为受温度制约.由于冰层内部是否发生热对流仍然没有定论 (Nimmo, 2004; Showman and Han, 2004), 本文假设能量以热传导的方式传递, “计算温度场”中将求解热传导方程.在实验数据约束下, 将冰视为非牛顿流体, “计算速度场”中将求解Stokes方程得到速度场.由速度场更新计算网格.判断当前 (t) 陨石坑深度h_t是否小于初始深度h₀的1/e, 如果是, 则停止计算; 如果否, 则进行下一步 (t+dt) 计算.下文将对计算细节分节描述.

图 1

Fig. 1

图 1 本文有限元热-力耦合计算流程图 Fig. 1 Calculation flow chart of FEM thermal-mechanical coupling

陨石撞击之后, 在黏性松弛作用下, 陨石坑形状将缓慢恢复 (Scott, 1967; Parmentier and Head, 1981).相对于岩石, 冰的黏度很小, 相应的松弛时间较短.这意味着, 相对于冰, 岩石能更好地保存初始陨石坑的形状.前人在研究冰体陨石坑时, 根据新近的月球陨石坑形状, 建立陨石坑初始几何模型 (Thomas and Schubert, 1986, 1987, 1988).依此思路, 本文陨石坑初始几何形状如图 2所示, 实线内为计算区域, h₀为陨石坑初始深度, D为陨石坑直径.若松弛时间相同, 地形波长越大, 则对应的冰层越厚 (Nimmo, 2004).欧罗巴陨石坑的直径最大约为40 km (Schenk, 2002).为了探求冰层厚度的最大值, 本文陨石坑直径取40 km.本文将针对不同厚度的冰层 (H=5km, 10km, 20 km) 分别计算陨石坑的松弛过程.

图 2

Fig. 2

图 2 几何模型及边界条件 实线内为计算区域; D=40 km; 不同模型冰层厚度H不同.温度边界条件:顶面温度Ts=95 K, 底面温度Tb=260 K, 侧面为绝热条件; 速度边界条件:顶面为自由面 (free surface), 侧面与底面为自由滑动面 (free slip). Fig. 2 Geometric model and boundary conditions

本文采用四边形网格对计算区域进行剖分.陨石坑附近的计算网格精度 < 0.1 km, 其余网格精度 < 0.5 km.

不考虑生热项的暂态热传导方程如式 (1) 所示:

(1)

其中c为比热, ρ为密度, T为绝对温度, t为时间, k为热导率.对于冰, c=2.05 (kJ·kg^-1·K^-1), ρ=916.5 (kg·m^-3), k=567/T(W·m^-1·K^-1)(Klinger, 1980).边界条件:冰层顶面温度T_s=95 K (Allu Peddinti and McNamara, 2015); 冰层底面温度T_b=260 K (Nimmo, 2004); 两侧为绝热 (如图 2).

在流体力学范畴内，地质材料 (包括冰) 可当作流体对待 (Karato, 2008).冰体的力学行为由方程 (2) -(4) 控制:

(2)

(3)

(4)

其中, μ_e为有效黏度 (Pa·s), v为速度 (m·s^-1), p为压力 (Pa), g=(1.3, 0) ^T为欧罗巴星的重力加速度 (m·s^-2), σ为偏应力张量 (Pa), 为应变速率张量 (s^-1).公式 (2) 是动量守恒方程, 即Stokes方程; 公式 (3) 是质量守恒方程, 本文假定冰为不可压缩流体; 公式 (4) 是本构方程.边界条件:顶面为自由面, 底面及侧面为自由滑动面 (如图 2).

本文视冰为非牛顿流体, 相应的差应力σ_d和应变速率的关系 (Goldsby and Kohlstedt, 2001) 如式 (5) :

(5)

其中A, n, p为流变常数, d为颗粒大小, Q为激活能, R为气体常数 (8.314 J·K^-1·mol^-1).由于缺少欧罗巴冰颗粒大小的实测数据, 本文根据地球冰体的数据推测.典型大陆冰颗粒大小为1~10 mm (Budd and Jacka, 1989).假设欧罗巴冰颗粒大小为5.5 mm.颗粒大小决定着冰体的流变性质.冰体单轴蠕变实验表明, 粒径5.5 mm对应超塑性流体 (superplastic flow)(Goldsby and Kohlstedt, 2001), 相关流变参数取值: A=3.9×10^-3 MPa^-1.8m^1.4s^-1, Q=49 kJ·mol^-1, n=1.8, p=1.4.

在数值模型中引用单轴实验数据时, 需要做转换, 具体公式如式 (6) 所示, 推导过程参见 (Gerya, 2010).实际计算中联立式 (2), (3), (4), (6) 求解速度场.

(6)

对Stokes方程不可压缩性的约束, 目前广泛采用的方法有两种, 即拉格朗日乘子法和罚函数法 (Bathe, 1996; Hughes, 2000).众多地球科学研究者采用罚函数法 (Dabrowski et al., 2008; Ismail-Zadeh and Tackley, 2010; Yang and Shi, 2015).本文拟采用此方法求解Stokes方程.关于罚函数方法的起源及可靠性问题的数学细节, 可参考 (Reddy, 1982; Cuvelier et al., 1986).

质量守恒方程 (3) 式的罚形式如 (7) 式所示, 其中λ为罚参数 (可理解为体积黏度).罚参数的引入使得材料具有一定的可压缩性.为了满足 (3) 式不可压缩的假设, 取足够大的罚参数, 这样就可以用 (7) 式代替 (3) 式进行有限元计算.一般地, 罚参数比剪切黏度大6~7个数量级 (Hughes et al., 1979; Donea and Huerta, 2003).将 (7) 式代入 (2) 式, 消除压力项, 便得到 (8) 式.这样, 原始的鞍点Stokes问题变为了消除压力项的Stokes方程 (Benzi et al., 2005).

(7)

(8)

目前, 有限单元法是求解偏微分方程最强大的数值方法之一 (Hughes, 2000; Dhatt et al., 2013).方程 (8) 对应的伽辽金有限单元公式如 (9) -(12) 式所示.其中，N为形函数矢量；B为形函数的空间偏导矩阵.通过求解 (9) -(12) 式, 我们将得到速度分布.

(9)

(10)

(11)

(12)

得到速度分布以后, 采用Uzawa迭代 (Benzi et al., 2005) 求解压力.该迭代法的基本算法如下：

(1) 选择任意的初始压力p⁰.

(2) 求解方程 (13).其中k为迭代数.

(3) 由 (14) 式校正压力.

(4) 重复第2步和第3步, 直到速度差满足一定的条件, 停止迭代.

(13)

(14)

陨石坑及附近地形随时间会发生缓慢变化.本文运用更新网格的办法示踪陨石坑表面形变.首先, 通过速度场更新表面计算节点的空间坐标.其次, 更新计算区域内部节点的空间坐标.我们曾成功地利用该方法追踪三维变形过程, 具体细节参见 (Yang and Shi, 2015).

欧罗巴冰层厚度制约着陨石坑松弛过程及松弛时间.本文针对不同的冰层厚度，计算了陨石坑的松弛过程，结果如图 3所示.黏性特征控制着力学行为.公式 (6) 中η₀(T) 可以反映冰层的黏性特征.η₀仅随温度变化而变化.温度场本身受传导机制控制, 所以比较稳定.基于上述, 我们考察不同厚度冰层模型初始时刻η₀的特征, 进而分析冰层厚度对松弛过程的制约.初始时刻陨石坑附近的黏度η₀如图 4所示.

图 3

Fig. 3

图 3 陨石坑深度随时间的变化图 h为陨石坑深度; t为时间; H为冰层厚度; h0为陨石坑初始深度; 30~80 Ma (阴影区) 为欧罗巴表面构造的平均年龄范围. Fig. 3 Temporal variation of impact crater depth

图 4

Fig. 4

图 4 初始时刻陨石坑附近的黏度η0 (a)(b)(c) 的冰层厚度分别为5 km, 10 km, 20 km. Fig. 4 Viscous coefficient nearby impact crater at starting time

图 3显示了不同冰层厚度下陨石坑深度h随时间的变化.图中可见: 1) 陨石坑深度随时间逐渐减小; 2) 陨石坑深度与时间的对数值呈近线性关系; 3) 冰层越厚, 陨石坑深度衰减越慢, 所需松弛时间越长 (5 km, 10 km, 20 km厚的冰层对应的松弛时间分别为0.85 Ma, 6.73 Ma, 49.81 Ma); 4) 欧罗巴表面构造的平均年龄介于30~80 Ma (Zahnle et al., 2003)(阴影区), 冰层厚度为20 km的模型对应的松弛时间落入此时间范围.这意味着, 如果冰层厚度小于20 km, 由于冰的黏性松弛, 欧罗巴表面存在直径为40 km陨石坑的可能性很低.鉴于欧罗巴陨石坑直径最大为40 km的事实 (Schenk, 2002), 本文认为冰层厚度大于20 km.值得注意的是, 我们和前人 (Thomas and Schubert, 1986, 1987, 1988; Nimmo, 2004) 的时间标准均为松弛时间, 如果改变此时间标准, 结论会有一定差异.

本文研究思路和Nimmo (2004)相近, 而研究结果不同:本文认为欧罗巴冰层厚度大于20 km, 而Nimmo (2004)认为小于10 km.结果有较大差异的原因可能是Nimmo (2004)存在疏漏.原文公式 (2) 中指数的分母缺少常数c, 更重要的是, 作者在引用单轴压缩实验数据时未作类似本文公式 (6) 的转换.相关转换公式的推导过程可参见Gerya (2010).

由公式 (6) 可知, η₀与exp (1/T) 呈线性关系, 那么log₁₀(η₀) 与1/T呈线性关系.模型中温度随深度近线性升高, 而图 4显示log₁₀(η₀) 随着深度近线性降低, 这与上述分析相符.另一方面, 在5 km深度内, 不同模型同一深度处log₁₀(η₀) 值不同:冰层较薄的模型log₁₀(η₀) 较小.这意味着, 对于大小相同的陨石坑, 冰层越薄则越容易发生变形, 相应的松弛时间越短, 这与图 3所示的结果相同.可见, 真正控制冰层流变特性的变量是温度.

欧罗巴表面存在许多小陨石坑 (直径小于40 km), 厚度大于20 km的冰层是否支持这些小陨石坑的存在呢?如果按照本文陨石坑形状推演, 小陨石坑的深度比大陨石坑的深度小.那么, 冰层厚度为20 km的情况下, 小陨石坑附近冰体的黏度将更大, 形状越难发生变形, 松弛时间将更长.换句话说, 厚度大于20 km的冰层支持小陨石坑的存在.

欧罗巴冰层下基岩面的温度可能对本文结果有一定影响.人们对于欧罗巴冰层之下基岩面的温度知之甚少, 作为一个端元模型, 本文跟随前人 (Nimmo, 2004) 取不会使冰体融化的温度 (260 K).另一个端元模型是基岩面的温度足以融化冰体.这意味着表面冰层与基岩之间将会出现水层.显然, 水层的存在将影响陨石坑的松弛过程.我们进一步的研究中将考虑水层的存在对陨石坑松弛过程的影响.目前, 有学者认为欧罗巴存在水层 (Carr et al., 1998; Kivelson et al., 2000; Spohn and Schubert, 2003).水层是广泛存在还是局部存在?此问题暂没有定论.鉴于上述, 本文认为实际的情况可能是部分冰层下有水层而部分冰层直接与基岩接触.相对于前述的两个端元模型, 此中间模型中陨石坑的松弛过程如何?我们进一步的研究中将深入分析.

本文基于有限单元法, 将冰视为非牛顿流体, 依据冰体实验数据, 修正了前人研究的不足, 对不同厚度冰体下陨石坑的松弛过程进行了动力学模拟, 进而估计了欧罗巴冰层的厚度.模拟结果显示: 1) 冰层越厚, 所需松弛时间越长; 2) 冰层越厚, 陨石坑附近的黏度越高, 这是松弛时间相对较长的直接原因.本文认为欧罗巴冰层的厚度大于20 km.值得注意的是, 作为端元模型, 本文模型中冰层与基岩直接接触.后续研究将进一步考虑其他模型.

两位匿名审稿人提出了建设性的意见, 编辑们付出了辛勤的劳动, 一并致以真挚的谢意！