ACTA Scientiarum Naturalium Universitatis Pekinensis
Numerical Study on Effects of Lateral Variations of Moon Crustal Thickness on Lunar Seismic Wave Propagation
CHEN Fei, WANG Yanbin†
Department of Geophysics, School of Earth and Space Sciences, Peking University, Beijing 100871; † Corresponding author, E-mail: ybwang@pku.edu.cn
Abstract A 2-D staggered grid pseudospectral and finite difference hybrid method is applied to perform numerical simulations of lunar seismic wave propagation in a laterally heterogeneous crustal Moon model. The relief of lunar crust and mantle interface is defined with different height, width and position based on the present lunar structural models. The effects of these factors on lunar records are discussed, respectively. How the lateral variations of Moon crustal thickness affects the strong wave coda of lunar signals is revealed. The results show that the height of the relief dominates the formation of strong wave coda, greatly affects the duration time and amplitude of the coda. The width of the relief has impact on the amplitude. The effect of position of the relief is complex, not only affecting the travel time of P and S waves, but also the duration time and amplitude of the coda. Key words lunar structure; lunar seismic wave propagation; numerical modeling; Moon crustal thickness variations
人类对月球的起源及演化充满兴趣, 这不仅因为月球是地球唯一的天然卫星, 更因为月球保存着早期的地质演化记录, 而这样的记录在地球上已经消失[1]。作为类地天体, 月球是极佳的研究样本。随着现代科学技术的发展, 人类对月球的了解已经从早期的平均半径、平均密度和质量等基本物理性
质发展到矿物组成、地形起伏以及重力场等更复杂
[2–3]和详细的信息 。然而, 对于研究月球的内部结构来说, 这些资料远远不够, 还需要获得月球内部的观测资料, 特别是月震波的观测资料。地震波是研究地球内部结构的有效方法之一, 同理, 月震波对研究月球内部结构具有非常重要的作用[4]。
国家自然科学基金(41374046)资助收稿日期: 20170223; 修回日期: 20170418; 网络出版日期: 20170510
迄今为止, 人类获得的月震波资料全部来自美国 Apollo登月计划在月球表面布设的月震仪。自1969 年 Apollo 11 号飞船首次在月球表面成功地布设月震仪以来, 月震仪在月球表面工作了 8 年时间,记录到累计超过 12000 次的月震[5]。在这些月震资料中, 最具科研价值的来自分别由 Apollo12, 14, 15和 16 号飞船布设的4台月震仪组成的观测台阵。台阵呈近似等边三角形, 每条边的长度约为 1200 km,其中 Apollo12 和 14 号飞船布设的月震仪大致位于三角形的同一个顶点, 间距仅 180 km。
Apollo 计划观测台阵记录的月震资料主要可以分为四大类: 天体撞击引发的月震、热成因月震、浅源月震和深源月震。与典型的地震记录上通常可以分辨出多个震相不同, 月震事件初至信号振幅很小, 事件持续时间长, 衰减缓慢[6], 后续震相被强烈的尾波掩盖。持久而强烈的尾波是月震信号不可忽视的典型特征。这种现象到底是什么原因导致的, 目前学界尚无定论。Blanchette-guertin 等[7]认为, 这种现象主要与月壳上部粗风化层的厚度有关。类似的现象也存在于区域地震记录中, 较多见于地震波在大陆地壳内部传播的过程, 一种可能的形成机制是地壳厚度的变化。
基于重力数据反演的研究结果[8]表明, 月壳厚度存在强烈的横向非均匀性, 特别是在撞击坑周围。与地壳厚度相比, 月壳厚度横向变化非常剧烈,在撞击坑内部, 最薄的月壳厚度接近 0 km, 月幔直接暴露于月表, 而月壳最厚的地方厚度超过 60 km。因此, 一个合理的推测是, 月壳厚度横向变化是月震波发生强烈散射与多次反射, 进而导致强烈尾波的重要因素之一。
对于地壳厚度变化对地震波的影响, 有过较多研究。Ruzaikin 等[9]讨论了地壳厚度横向变化对 Lg波传播的影响。Zhang 等[10]通过波形模拟和观测资料的对比, 讨论了地形起伏对于 Lg 波和 Sn 波传播的影响。Furumura 等[11–12]采用有限差分数值模拟,研究俯冲带地壳从陆壳向洋壳过渡对 Lg 传播的影响。Liu 等[13]采用谱元法, 模拟地壳厚度变化对瑞雷面波到时的影响。此外, Larmat 等[14]采用模式算法与谱元法相结合的方法, 模拟火星地形, 研究火星地壳与地幔分界面起伏对瑞雷面波传播的影响。总结起来, 无论是针对地球或者火星, 在研究地壳厚度变化对地震信号的影响时, 都是关注某一具体震相, 比如 Lg波、瑞雷面波等。
[15] Jiang 等 采用有限差分和伪谱混合方法, 模拟全月球模型中 SH 波的传播, 讨论月壳厚度阶梯状横向变化对月震信号的影响。Blanchette-guertin等[7]采用改进的声子法, 模拟存在散射的情况下阶梯状盆地月壤层厚度变化对月震波传播的影响。声子法的原理基于射线理论, 作为一种间接方法, 声子法虽然计算效率高, 但包含的波场信息不如直接求解波动方程那样完整和准确。Jiang 等[15]的月震模拟虽然基于波动方程直接求解, 但只考虑了月壳厚度变化对月震记录的影响, 不具有一般性。
本文选择 Jiang 等[15]采用的基于交错网格的伪谱与有限差分混合的方法, 根据现有月球结构研究结果, 设计更接近实际的撞击导致月壳厚度横向变化的起伏模型, 模拟不同起伏高度、宽度和位置对月震波传播的影响。
1研究方法与模型1.1交错网格伪谱与有限差分混合方法
考虑P波和SV波在柱坐标(r, , z)下的传播。如果假设 Z 方向上的所有变量都不变, 那么可以得到 P-SV 波在二维柱坐标下一阶应力速度形式的波动方程:
其中, vr 和 v分别表示质点在径向和横向的运动速度, 表示介质密度, fr 和 f 分别表示径向和横向上的体力, , 和 表示应力分量。对于各向同性弹性介质, 应力与速度的关系为
其中,为了求解方程(1)和和表示拉梅常数。(2), 采用如图 1所示的交错网格, 在r和 方向上对模型进行离散化。速度和应力分量在相隔半个网格节点的不同位置进行离散化。在 方向上, 对于给定的半径, 采用相同的网格间距 进行离散; 在 r方向上, 采用相同的网
格间距r进行离散。采用基于交错网格的伪谱与有限差分混合方法, 计算方程(1)和(2)中的空间偏导数。在横向上,采用快速傅里叶变换计算对 的偏导数[17]:
方程(1)和(2)中, 计算对时间的偏导数时, 采用与Furumura 等[19]相同的有限差分形式。
我们选取通过月球大圆的深度为 1000 km 在水平方向上延伸 120°的扇形剖面为模拟月震波传播的区域(图 2)。在扇形区域中需要处理 4 个边界条件: 对于月球表面, 采用牵引力为零的自由边界条件; 为了减少人为反射的影响, 对于左右两侧和底边界, 采用 Cerjan 等[20]给出的 20 个节点吸收边界条件。方程(1)中的震源体力 fr 和 f 采用 Wang 等[16]给出的二维柱坐标下的线源形式: 其中, r0 和 为震源的中心位置, Mrr , Mr 和 M 0 分别为震源矩张量分量。式(5)和(6)中的 函数采用 Herrmann 伪 函数近似。Jiang 等[15]和 Wang 等[16]分别用该二维模型计算了全月球模型中 P-SV 和 SH 波的传播, 结果可以与实际月震观测波形进行对比, 说明了该二维模型的有效性。在本文的横向非均匀月壳模型中, 考虑到横向非均匀的尺度远大于计算的地震波波长,因此, 二维模型仍然可以保证计算的有效性。1.2 月球内部结构模型[21] [22]基于Khan等 、Lognonné等 和 Gagnepainbeyneix [23]等 针对月球一维速度结构的研究结果, [24] Weber 等 2011 年采用阵列处理的方法, 对现有的 Apllo 数据重新进行分析, 发现来自月球内核的反射波和转换波, 并提出月球外核为液态的观点。[24]本文数值模拟中, 采用 Weber 等 给出的月球分
层速度和密度结构模型作为背景模型(图 3(a))。
[25] Garcia 等 2011年对相关月震重新定位, 发现来自月核反射的 SH 波。他们通过对月震数据和重力数据的联合反演, 得到非常初步的月球参考模型(VPREMOON-VERY Preliminary Reference Moon Model)。本文数值模拟中, 采用 Garcia 等[25]给出的滞弹性衰减Q值作为参考值(图 3(b))。
与地球相比, 月球表面遭受过强烈而持续的天体撞击, 导致月球表面形成大小不一的撞击坑, 半径从数千米至数百千米。一般将半径超过 200 km的撞击坑称为盆地。这些撞击坑的存在, 正是导致壳均衡假说解释陨石坑下剧烈的月壳‒月幔分界面月壳厚度横向剧烈变化的最重要原因。有人曾用地
[8]起伏。Wieczorek 等 基于重力数据反演的研究结果表明, 在撞击坑内部, 最薄的月壳厚度接近 0 km, 月幔直接暴露于月表, 而月壳最厚的地方厚度可超过 60 km。基于 Hikida等[26]给出的月球盆地周计了正弦形状的月壳‒月幔分界面起伏,围月壳厚度剖面图, 在背景模型的基础上, 我们设
定义了起伏的 3 个要素: 高度、宽度和震中距(起伏中心与震中的距离)(图 2)。这是比简单的阶梯状月壳厚度变化更接近真实的月壳模型。
2 数值模拟研究
本节通过数值计算, 模拟月壳伏三要素对浅震和深震月震波传播的影响。月球内 部的浅震震源深度一般为 50~220 km; 深震具有丛集性, 震源深度通常为 700~1200 km。在本节的模拟中, 设计浅震深度为 100 km, 深震为 900 km。浅震和深震的震源机制均采用二维的双力偶线源, 取Mrr 1.0, M 1.0, Mr M 0 , 相当于 45°倾
r角的走滑断层。
综合考虑计算机性能和计算效率, 对图 2 所示的扇形研究区域进行如下的网格剖分。将 方向离散为 2048 个节点, r方向离散为 1000 个节点。如此, 方向的网格间距最大为 1.78 km, 位于模型表面, 最小为 0.75 km, 位于模型底部。r方向的网格间距为 1 km。根据 Wang 等[16]和 Furumura 等[19]的研究, 考虑最大网格的分辨能力, 模拟中的最小波长 4.53 km。考虑到月球模型中的最小波速
min Vmin=1.8 km/s, 可以得出所能模拟的最高频率为f Vmin / min 0.4 Hz。因此, 时间震源函数采用max宽度为 2.5 s 的 Herrmann 伪 函数。为了保证计算的稳定性, 时间步长 t 应当满足稳定性条件t / Vmax , 其中 为常系数 0.26, 故模拟中
min ‒月幔分界面的研究成采用 t 0.02 s。本文根据前人关于月壳果, 设计的分界面起伏宽度分别为 180, 360 和 720 km, 尺度远远大于模拟的月震波波长。所以, 从“Z方向上所有变量都不变”推出二维模型的假设条件,在全球尺度上依然成立。
2.1 起伏高度对月震波传播的影响
‒月幔分界面起伏宽度为
将震源位置放在 120°的扇形模拟区域内 20°处(图 2), 固定月壳 360 km,起伏中心在 60°处, 设计 5 个不同高度的起伏模型,分别为 0, 10, 19.5, 30 和 39 km (见表 1), 研究起伏高度对浅震(100 km)和深震(900 km)月震波传播的影响。
从图 4 看出, 当起伏高度逐渐增加时, 会引起越来越强烈的多次反射波和各种转换波, 它们的相干叠加会形成次生面波。这些波的共同作用造成强烈而持久的尾波, 在水平分量上表现得尤为明显,
其传播速度明显低于直达波, 其成因类似地球上沉积盆地边缘产生的多重反射、转换和次生面波。当从平层月壳结构变为有 10 km 的高度起伏时, 月震信号明显受到影响。可以看出, 起伏中心(震中距40°)之后的月震波在 ss 震相之后的信号振幅有所增加, 并且持续时间增长, 但 P, pp, S 和 ss 等震相未受影响。当起伏高度增至 39 km 时, 这种影响更加明显, 尾波振幅和持续时间急剧增加, 起伏之前的月震波也受到影响。近震源一侧的起伏会阻挡月震波在月壳中的传播, 将信号反射回去, 在震中距−20°~40°范围内可以清晰地辨别此反射波, 其波速与震中距 40°~100°的尾波相同, 但是振幅较弱。
图 5 显示震中距 60°处月震波在不同起伏高度下的水平分量。可见, 不论是浅震还是深震, 由于起伏的影响, 在起伏之后的月震信号均会产生尾波,并表现出明显的低频面波特征。对于浅震, 大约700~900 s 处的月震波振幅明显增大。当起伏高度增加到39 km 时, 900~1500 s 之间出现强烈而持久的大振幅尾波。对于深震, 月震波的整体振幅比浅震信号小, 但仍可清楚地看到 450~750 s 之间, 月震波随着起伏高度的增加振幅逐渐增加, 当起伏高 度达到 39 km 时, 750~1500s 之间出现强烈而持久的尾波。P, pp, S 和 ss等震相的到时振幅均未受影响。
从上述讨论可知, 起伏高度主要影响月震波尾波的振幅和持续时间。
2.2 起伏宽度对月震波传播的影响
固定月壳‒月幔分界面起伏高度为将震源位置放在 120°的扇形模拟区域内 20°处(图 2), 39 km, 起伏中心在 60°处, 设计 3 个不同宽度的起伏模型,其宽度分别为180, 360和 720 km (见表 1), 研究起伏宽度对浅震(100 km)和深震(900 km)月震波传播的影响。
从图 6 可以看出, 当起伏宽度逐渐增加时, 起伏中心(震中距 40°)之后的月震波尾波持续时间未见明显变化, 但振幅有所增加, 而 P, pp, S, ss 等震相的到时振幅均未受影响。
图 7 显示震中距 60°处, 月震波在不同起伏宽度下的水平分量。可见, 不论是浅震还是深震, 月震波的尾波持续时间基本上保持不变, 但振幅有所增加。对于浅震, 当起伏宽度增加到 720 km 时, 600~1200 s 之间的月震波振幅明显增加, 而深震信号振幅的增加主要体现在 1200~1500 s 之间。随着
起伏宽度的逐渐增加, 月壳厚度会逐渐变薄, 月壳内反射波的反射和转换次数随之增加, 使得多重反射波与转换波更加密集地叠加, 从而导致尾波的振幅增加。
从上述讨论可知, 起伏宽度主要影响月震波尾 波的振幅。固定月壳‒月幔分界面起伏高度为2.3 起伏位置对月震波传播的影响39 km, 起伏宽度为 360 km, 从 20°, 40°到60°移动震源位置(相当于起伏中心与震源之间的距离从 40°减小到 20°,
0°, 见表 1), 研究起伏位置对浅震(100 km)和深震(900 km)月震波传播的影响。
从图 8 可以看出, 起伏位置对于月震波传播的影响很大, 但在浅震和深震月震波传播中, 均未见明显的规律。
图 9 显示震中距 60°处, 月震波在不同起伏位 置下的水平分量。当起伏从震中距 40°减小到 20°再到 0°时, 月震波变化明显。P, S 等震相到时提前,月震波尾波的持续时间和振幅变化剧烈。对于浅震, 起伏中心位于震中距 40°时, 月震波的尾波持续时间和振幅均大于起伏中心位于 20°下的月震信号。但是, 当起伏中心位于震中距 0°时(即起伏中
心位于震源正上方), 会产生强烈而持久的大振幅月震尾波。对于深震, 起伏中心位于震中距 40°时月震波的尾波持续时间和振幅却小于起伏中心位于震中距 20°时的月震信号。当起伏中心位于震中距0°时, 月震尾波的持续时间增长, 但振幅减小。起伏的位置主要影响月震波的入射角度。月震波穿过起伏界面时, 是从波密介质入射到波疏介质, 此时入射角的变化对透射波的能量影响不大。但是, 当透射波在月壳中多次反射与转换时, 则从月壳透射至月幔的能量与入射角度的关系密切, 但并非简单的线性变化, 由此导致复杂变化的月震尾波。
综上所述, 起伏位置不仅影响 P, S等震相的到时, 而且影响月震波尾波的持续时间和振幅。
3 结论和讨论
‒月幔分界面起伏模型。模本文采用基于交错网格的二维伪谱和有限差分混合方法, 设计了月壳拟了浅震和深震月震波在不同起伏高度、宽度和位
月壳‒月幔分界面的起伏高度在月壳厚度横置下的传播。通过数值模拟, 得到以下结论。
1)向变化对月震波传播影响中起主导作用。起伏高度越大, 表示月壳厚度横向变化越剧烈, 会造成月震波强烈的多次反射和转换, 多次反射波和转换波的
月壳‒月幔分界面起伏宽度对月震信号传播相干叠加导致强烈的月震尾波。
2)也有影响, 在起伏高度一定的情况下, 越宽的起伏
月壳‒月幔分界面起伏中心相对于震源所在会造成越强的多次反射, 导致越强的月震尾波。
3)位置对月震波传播影响很大, 但难以预测。对浅震和深震信号的影响并不一致, 这与月震波相对于起伏界面的入射角有密切关系。
由于计算性能限制, 本文模拟的月震信号最高频率仅为 0.4 Hz, 而实际月震信号频率可高达 3 Hz以上。更高频的数值模拟将有助于更好地探究月壳厚度横向变化对月震波传播的影响。同时, 本文未讨论月表地形起伏, 然而实际月壳厚度横向变化的成因包含地形起伏。此外, 月壳三维横向非均匀结构也会影响月震波的散射。在后续研究中, 若能将地形因素以及三维非均匀结构考虑进去, 将进一步提高数值模拟结果的可靠性。
参考文献
[1] Konopliv A S, Binder A B, Hood L L, et al. Improved
1476‒1480 gravity field of the Moon from lunar prospector. Science, 1998, 281: [2] Lawrence D J, Feldman W C, Barraclough B L, et al. Thorium abundances on the lunar surface. Journal of 20307‒20331 Geophysical Research Atmospheres, 2000, 105(E8): [3] Feldman W C, Gasnault O, Maurice S, et al. Global distribution of lunar composition: new results from
5-1‒5-14 Lunar Prospector. Journal of Geophysical Research Planets, 2002, 107:
15‒24 [4] 姜明明, 艾印双. 月震与月球内部结构. 地球化学, 2010, 39(1): [5] Nakamura Y. Farside deep Moonquakes and deep interior
251‒268 of the Moon. Journal of Geophysical Research, 2005, 110(E1): [6] Dainty A M, Toksöz M N, Anderson K R, et al. Seismic scattering and shallow structure of the Moon in
11‒29 Oceanus Procellarum. Earth, Moon, and Planets, 1974, 9(1): [7] Blanchette-guertin J, Johnson C L, Lawrence J F. Effects of lateral variations in megaregolith thickness
171‒178 on predicted lunar seismic signals. Geophysical Research Letters, 2015, 42: [8] Wieczorek M A, Neumann G A, Nimmo F, et al. The
671‒675 crust of the Moon as seen by GRAIL. Science, 2013, 339: [9] Ruzaikin A I, Nersesov I L, Khalturin V I, et al. Propagation of Lg and lateral variations in crustal
307‒316 structure in Asia. Journal of Geophysical Research, 1977, 82: [10] Zhang T, Lay T. Effects of crustal structure under the Barents and Kara Seas on short-period regional wave propagation for Novaya Zemlya explosions: empirical
1132‒1147 relations. Bulletin of the Seismological Society of America, 1994, 84(4): [11] Furumura T, Kennett B L N. Subduction zone guided waves and the heterogeneity structure of the subducted plate: intensity anomalies in northern Japan.
1‒11 Journal of Geophysical Research Atmospheres, 2005, 110(B10): [12] Furumura T, Hong T K, Kennett B L. Lg, wave propagation in the area around Japan: observations and
1‒20 simulations. Progress in Earth and Planetary Science, 2014, 1(1): [13] Liu K, Zhou Y. Effects of crustal thickness variations on surface wave phase delays. Geophysical Journal
773‒792 International, 2013, 192(2): [14] Larmat C, Montagner J P, Capdeville Y, et al. Numerical assessment of the effects of topography and crustal thickness on martian seismograms using a
78‒89 coupled modal solution-spectral element method. Icarus, 2008, 196(1): [15] Jiang X, Wang Y, Qin Y, et al. Global Sh-wave propagation in a 2D whole Moon model using the
163‒174 parallel hybrid PSM/FDM method. Earthquake Science, 2015, 28(3): [16] Wang Y, Takenaka H, Jiang X, et al. Modelling twodimensional global seismic wave propagation in a
1271‒1287 laterally heterogeneous whole-moon model. Geophysical Journal International, 2013, 192(3): [17] Witte D C, Richards P G. The pseudospectral method
1‒18 for simulating wave propagation. Computational Acoustics, 1990, 3:
1425‒1436 [18] Levander A R. Fourth-order finite-difference P-SV seismograms. Geophysics, 1988, 53(11): [19] Furumura T, Kennett B L N, Furumura M. Seismic wavefield calculation for laterally heterogeneous
845‒ whole earth models using the pseudospectral method. Geophysical Journal International, 1998, 135(3): 860 [20] Cerjan C, Kosloff D, Kosloff R, et al. A nonreflecting boundary condition for discrete acoustic and elastic
705‒708 wave calculations. Geophysics, 1985, 50(4): [21] Khan A, Mosegaard K, Rasmussen K L. A new seismic velocity model for the Moon from a Monte Carlo inversion of the Apollo lunar seismic data. Geophysical Research Letters, 2000, 27(11): 1591– 1594 [22] Lognonné P, Gagnepain-beyneix J, Chenet H. A new seismic model of the Moon: implications for structure, thermal evolution and formation of the Moon. 27‒44 Earth & Planetary Science Letters, 2003, 211(1/2): [23] Gagnepain-beyneix J, Lognonné P, Chenet H, et al. A seismic model of the lunar mantle and constraints on
140‒166 temperature and mineralogy. Physics of the Earth & Planetary Interiors, 2006, 159(3/4):
309‒312 [24] Weber R C, Lognonné P. Seismic detection of the lunar core. Science, 2011, 331: [25] Garcia R F, Gagnepain-beyneix J, Chevrot S, et al.
96‒ Very preliminary reference Moon model. Physics of the Earth & Planetary Interiors, 2011, 188(1/2): 113 [26] Hikida H, Wieczorek M A. Crustal thickness of the
150‒ Moon: new constraints from gravity inversions using polyhedral shape models. Icarus, 2007, 192(1): 166