2. 中国科学院大学, 北京 100049
3. 中国地震台网中心, 北京 100045
2. University of Chinese Academy of Sciences, Beijing 100049, China
3. China Earthquake Networks Center, Beijing 100045, China
地球自由振荡是研究地球内部物质结构和特征的重要信息来源之一.18世纪初,研究者就已经开始了对地球自由振荡的理论研究.1829年,Poisson探讨了完全弹性固体球的自由振动问题(Shearer, 2009).1863年Kelvin发展了弹性球体的应变理论并将其应用于固体潮,开创了对地球自由振荡的研究(Lapwood and Usami, 1981).1881年,Lamb解决了均匀弹性球在无重力无自转情况下的振动问题,并提出了自由振荡的基本类型(Lapwood and Usami, 1981).随着对地球自由振荡的深入研究,同行们一步步克服了真实地球构造所带来的巨大困难,对理论地球模型加入了更多的物理形态.1911年,Love基于重力作用讨论了可压缩球体的自由振荡现象(Bolt, 2012).Dahlen(1968)利用扰动理论,分析了地球自转及椭率作用导致的地球自由振荡谱分裂;基于这一成果,Rogister(2003) 进一步从基本模式理论出发,得出导致自由振荡谱分裂的主要原因是自转的一阶效应和科里奥利力的二阶效应.
由于自由振荡的波长很长且尺度很大,影响地球自由振荡频谱特征的因素有很多,比如地球结构的横向及径向不均匀性、自转、椭率、重力以及震源机制等因素,因此研究地球自由振荡的频谱特征,可以得到关于地球内部甚至是深部结构的详细信息,这对于研究地球内部非均匀结构以及地球自身物理状态等具有重要意义.
核幔边界D″区的横向非均匀性是探究地球深部结构与性质的重点研究内容之一.现有研究表明,D″区位于地幔底部、核幔边界以上250 km左右的范围,存在强烈的横向不均匀性,并会对地球的自由振荡产生一定影响;D″区下部存在一类极为异常的横向非均匀结构——超低速异常带,其厚度一般在5~60 km之间,横向尺度约200 km,剪切波速异常最高可达-30%,压缩波速异常可达-10%,它对D″区横向非均匀的演化起到至关重要的作用(倪四道,2007).
本文基于PREM地球模型,在核幔边界引入具有横向非均匀结构的D″区与超低速异常带,构造了新的地球模型.由此进行全球自由振荡数值模拟实验,即求解地球自由振荡波动方程,并结合大规模并行计算机模拟地震波的传播.数值模拟结果的对比分析表明,核幔边界横向非均匀结构以及不同地球模型的选取对自由振荡频谱有明显影响.
1 地球自由振荡理论与波动方程求解Alterman等(1959)把二阶偏微分方程本征值问题转化为一阶常微分方程组,由此拉开了计算机数值模拟的序幕.方俊(1990a,b)讨论了地球内部弹性参数和密度参数的变化对地球自由振荡本征周期的影响.Tsuboi等(1985)探究了地球自由振荡谱峰与横向非均匀性之间的可能关系,得出谱峰的中心频率随时间而改变的结论.严珍珍等(2008)采用谱元法,结合并行计算数值模拟特大地震激发的弹性波在地球内部的传播过程.
地球自由振荡可分为两类基本振型,一类是环型振荡,各质点只在以地心为球心的同心球面上振动,位移无径向分量,地球介质只产生剪切形变,没有体积变化,由nTlm表示;另一类是球型振荡,其质点位移既有径向分量,也有水平分量,由nSlm表示(Lapwood and Usami, 1981).
地球自由振荡的波动方程求解,类似于地震波的波动方程求解,本质上是求满足适当初边值条件的地球振动的常系数偏微分方程组的数值解(陈世仲等,2016).
在地壳与地幔结构中,考虑地球自转,忽略重力,引入科里奥利力的三维弹性波动方程为(Komatitsch and Tromp, 2002b):
(1) |
其中,Ω为地球自转角频率矢量,σ为应力矢量,f为体力项.
地表可看作自由表面,在自由表面上,边界条件应满足n·σ=0,其中,n是自由表面的法向量.在核幔边界(CMB)位移的法向量n·n应连续,且地幔底部的应力n·σ应与外核顶部的应力-pn一致,p是外核流体压力.用M表示地幔和地壳区域,则(1) 式的弱积分形式为
(2) |
在液态外核中,利用科林近似,同样忽略重力,波动方程可以写为
(3) |
其中,κ是流体的体积模量.引入一个标量势χ,则流体压力p=-κ∂·u也可以写为
(4) |
定义一个矢量函数φ,将位移场u表述为标量势函数χ和矢量函数φ的和形式,公式为
(5) |
最终,得到只考虑自转的波动方程为
(6) |
将(4) 与(5) 式代入(3) 式,令
(7) |
其中,系数A与B可由下两式表达,公式为
(8) |
矢量函数φ没有z分量,
(9) |
在内核中,视内核为弹性的固体结构,其弱形式与(2) 式相似,但内外核边界的应力、法项位移与速度应保持连续性.故其弱积分形式为
(10) |
在得到分层结构的理论地震波波动方程的弱积分表达式后,即可对其进行空间域离散和时间域离散.对于空间域离散,本文采用被广泛认同的谱元法(Boyd, 2000),根据谱元法的原理,需求解的方程应满足空间离散后的微分方程,公式为
(11) |
其中,U是所有网格点的位移矢量,M为对角质量矩阵,W是科里奥利力矩阵,K代表刚度矩阵,B为核幔边界和内外核边界的耦合项.
将波动方程运用谱元法进行空间域离散后,在时间域使用Newmark算法(Newmark, 1959)进行离散,即可求解.
2 数值模拟分析本文的数值模拟研究采用2012年4月11日16时38分(北京时间)发生在印度尼西亚苏门答腊西北海域(92.98°E,2.35°N)的8.6级强震参数,并以其地震记录、频谱及功率谱密度为基本参照,基于PREM地球模型分别引入核幔边界具有横向非均匀结构的D″区以及超低速异常带,构建两个新的横向非均匀结构地球模型,使用Computational Infrastructure for Geodynamics (CIG)提供的谱元程序包SPECFEM3D_GLOBE(Komatitsch and Tromp, 2002a, b),在“天河2号”高性能计算平台上利用600个进程进行并行计算.采用谱元法模拟地震波在全球传播500 min,根据地球自由振荡甚低频的特性,采用低通滤波0.0039 Hz以尽可能压制其他相关“干扰”因素.本文所有地球模型皆考虑地球自转.
PREM地球模型引入D″区的实现,为采用Kustowski(2007)提出的S362ANI三维地球模型核幔边界D″区的剪切波速度模型数据.图 1a所示为2800km深度处具有横向非均匀结构的D″区.根据Thorne等(2007)提出的南太平洋底部的大面积超低速异常带,我们假设这片区域剪切波速异常为-30%,如图 1b所示,深度为2890 km,图中红色区域为超低速异常带.
限于篇幅,在全球台站中仅选用PALK.Ⅱ(80.70°E,7.27°N)台站与MBWA.IU(119.73°E,21.16°S)台站展示数值模拟分析说明,并对数值模拟后的结果分别进行频谱分析与功率谱密度分析.
2.1 频谱分析分别将引入核幔边界D″区的PREM地球模型和引入超低速异常带的PREM地球模型的数值模拟数据进行滤波与归一化处理,得到时间域合成地震图.为了突出长时程模拟结果对比,选取200~500 min进行展示.如图 2所示,红色曲线为采用PREM地球模型,蓝色曲线为采用PREM地球模型引入D″区,绿色曲线为采用PREM地球模型引入超低速异常带时所得到的长时程理论地震图.利用chip-z变换对时间域合成地震图进行频谱转换,得到频率域以简正模为参考的频谱图,见图 3与图 4,为了突显模拟结果,截取0.0027~0.0039 Hz带宽进行对比.
从图 2可看出,虽然具有横向非均匀结构的D″区速度扰动幅度非常小,超低速异常带的尺度对于全球尺度来说亦不算大,但可看出二者在时间域上波形的明显变化,由此可知不论是全球尺度的横向非均匀结构,亦或是区域尺度的横向非均匀结构,皆对地球自由振荡有一定影响.
从图 3可看出,对于图中大多数振型,超低速异常带与D″区对地球自由振荡的影响均较明显;对于0S24、0S25、0S26和0S27振型,地球自由振荡对超低速异常带的敏感度明显高于具有全球尺度的D″区;且除0S27、0S28、0S29和0S30振型外,核幔边界的这两类横向非均匀对球型振荡基频振型的特征频率几无影响,仅对幅值有不同程度的影响.类似地,如图 4所示,对于0T21、0T24和0T25振型,地球自由振荡对超低速异常带的敏感度高于D″区;核幔边界地幔侧的横向非均匀对环型振荡基频振型的幅值有不同程度的影响,而对其特征频率均无影响.因此,基频振型中0S27、0S28、0S29和0S30振型适合于利用地球自由振荡谱线分裂现象分析核幔边界地幔侧横向非均匀结构的研究,而其余基频振型则不适合此类研究.
就频谱特征而言,在上述两类横向非均匀结构分别作用下,模型中地球自由振荡数值频谱曲线的变化形态相近,幅值略有差异.
2.2 功率谱密度分析为了更直观、定量地评估不同尺度横向非均匀对地球自由振荡的影响,采用功率谱密度估计的方法处理分别引入核幔边界D″区横向非均匀的PREM地球模型和引入超低速异常带的PREM地球模型的数值模拟数据,并对幅值进行定量分析,如图 5与图 6所示.图中红色曲线为采用PREM地球模型,蓝色曲线为采用PREM地球模型引入D″区,绿色曲线为采用PREM地球模型引入超低速异常带所对应的功率谱密度曲线图.这里只选取PALK.Ⅱ台站进行详细说明,MBWA.IU台站模拟结果与其大致相同.
由图 5可以看出,虽然三组数据的功率谱密度峰值与所对应的简正模并无较大差别,但横向非均匀结构对于自由振荡某些振型能量的影响较为明显.对此,我们进一步将幅值进行定量分析,如表 1所示.
表 1清楚地表明,PREM地球模型引入横向非均匀结构后振型间能量漂移明显,引入超低速异常带的振型能量耗散(功率谱密度幅值负偏差)比引入D″区的耗散更明显.由此可知,地球自由振荡对超低速异常带的敏感性高于对D″区.
由图 6与表 2可以看出,PREM地球模型引入横向非均匀结构后,幅值在T分量上并没有Z分量变化大.由此可见,就横向非均匀结构的影响而言,地球自由振荡的Z分量比T分量更敏感.
由上述功率谱密度分析可知,本文所采用的两个不同尺度的、核幔边界附近的横向非均匀结构对模型中地球自由振荡的特征频率的影响不明显,而对部分振型能量的衰减及各振型间的能量漂移影响较为明显.因此,类似于频谱分析的结论,除部分基频球型振型外,大部分基频振型都不适用于利用地球自由振荡谱线分裂现象分析核幔边界地幔侧横向非均匀结构的研究.另外,超低速异常带对模型中地球自由振荡各振型能量变化的影响明显大于横向非均匀D″区.
为了探究地球模型对模拟结果的影响,本文引入了PALK.Ⅱ台站的观测地震数据,在三维S362ANI地球模型(Kustowski,2007)的基础上,通过上述方法引入了超低速异常带,并与引入超低速异常带的PREM地球模型作比较.结果如图 7与图 8所示,图中红色曲线为采用PREM/S362ANI地球模型引入超低速异常带,蓝色曲线为采用观测数据所对应的功率谱密度曲线.
图 7与图 8可看出,采用包含横向非均匀结构的S362ANI三维地球模型自由振荡各振型的功率谱密度明显比PREM地球模型的更接近实际观测数据,这就意味着地球内部横向非均匀结构对地球自由振荡的影响是明显的、确定无疑的,故地球模型的选取对地球自由振荡的模拟至关重要.尽管如此,与实际观测数据相比,采用两种模型所得功率谱密度曲线均有明显偏差,这是因为两种模型皆未考虑地球介质衰减、重力、固体潮及地球扁率等影响因素,模拟结果与观测数据间必然存在明显差异.
3 结论与讨论 3.1 结论与讨论本文对2012年4月11日发生在印度尼西亚苏门答腊西北海域(92.98°E,2.35°N)的8.6级强震进行地球自由振荡波动方程数值模拟分析,在PREM地球模型的基础上分别引入了地球自转,并引入核幔边界具有横向非均匀结构的D″区与超低速异常带,采用了频谱分析与功率谱密度估计的方法进行分析.频谱分析与功率谱密度估计分析的结果表明,不同尺度的核幔边界横向非均匀结构对地球自由振荡均有一定影响.对于球型振荡,0S24、0S25、0S26和0S27振型对超低速异常带的敏感度明显高于D″区;对于环型振荡,0T21、0T24和0T25振型对超低速异常带的敏感度亦明显高于具有横向非均匀结构的D″区.这是由于D″区横向非均匀虽为全球尺度,但其扰动幅度非常小,约为2%~3%,远逊于超低速异常带的-30%扰动幅度.上述频谱分析与功率谱密度估计分析还表明,除0S27、0S28、0S29和0S30振型外,核幔边界的这两类横向非均匀对自由振荡大部分基频振型的特征频率几无影响,仅对幅值有不同程度的影响.这就意味着,基频振型中仅0S27、0S28、0S29和0S30振型适合于利用地球自由振荡谱线分裂现象分析核幔边界地幔侧横向非均匀结构的研究,而其余基频振型则不适合此类研究.
3.2 结论与讨论由于我们在模拟实验中所用模型的背景结构为一维模型,与实际地球结构相差较大,并且在进行波动方程空间离散时剖分的精细程度依然较粗糙.另外,在时间离散时采用的Newmark算法为非保结构算法,加之本文所采用之模型未考虑介质衰减、重力、固体潮及地球扁率等因素,引入横向非均匀结构的D″区和超低速异常带的PREM模型亦较为简化.后续研究将在三维模型的基础上实现空间网格剖分局部加密以及引入本课题组提出的保结构算法.由于所使用SPECFEM3D GLOB软件的剖分规则限制,本文没有实现对核幔边界地表起伏的模拟分析.
3.3 结论与讨论需要说明的是,本文所涉研究旨在评估核幔边界附近不同尺度(全球尺度与大区域尺度)的横向非均匀结构对地球自由振荡的影响,为后续相关研究提供参考和相关的反演约束;而非基于精确地球模型实现地球自由振荡的数值模拟.
致谢 感谢国家超级计算广州中心支持.[] | Alterman Z, Jarosch H, Pekeris C L. 1959. Oscillations of the earth[J]. Proc. R. Soc. Lond. A Math. Phys. Sci., 252(1268): 80–95. DOI:10.1098/rspa.1959.0138 |
[] | Boyd J P. 2000. Chebyshev and Fourier Spectral Methods[M]. 2nd ed. Mineola, New York: DOVER Publications. |
[] | Chen S Z, Li X F, Wang W S. 2016. Structure-preserving numerical simulation for Earth's free oscillations[J]. Chinese Journal of Geophysics (in Chinese), 59(5): 1685–1695. DOI:10.6038/cjg20160513 |
[] | Dahlen F A. 1968. The normal modes of a rotating, elliptical earth[J]. Geophysical Journal International, 16(4): 329–367. DOI:10.1111/gji.1968.16.issue-4 |
[] | Fang J. 1991a. The parameter kernels for the linear inversion of the free oscillation of the earth (Ⅰ)[J]. Science in China Series B, 34(1): 78–85. |
[] | Fang J. 1991b. The parameter kernels for the linear inversion of the free oscillation of the earth (Ⅱ)[J]. Science in China Series B (in Chinese), 34(1): 86–94. |
[] | Komatitsch D, Tromp J. 2002a. Spectral-element simulations of global seismic wave propagation-Ⅰ. Validation[J]. Geophysical Journal International, 149(2): 390–412. DOI:10.1046/j.1365-246X.2002.01653.x |
[] | Komatitsch D, Tromp J. 2002b. Spectral-element simulations of global seismic wave propagation-Ⅱ. Three-dimensional models, oceans, rotation and self-gravitation[J]. Geophysical Journal International, 150(1): 303–318. DOI:10.1046/j.1365-246X.2002.01716.x |
[] | Kustowski B. 2007. Modeling the anisotropic shear wave velocity structure in the earth's mantle on global and regional scales[Ph.D. thesis]. Cambridge:Harvard University. |
[] | Lapwood ER, Usami T. 1981. Free Oscillations of the Earth[M]. London: Cambridge University Press. |
[] | Newmark N M. 1959. A method of computation for structural dynamics[J]. Journal of the Engineering Mechanics Division, 85(3): 67–94. |
[] | Ni S D. 2007. Review on studies of ultra low velocity zone[J]. Journal of University of Science and Technology of China (in Chinese), 37(8): 820–829. |
[] | Shearer P M. 2009. Introduction to Seismology[M]. 2nd ed. New York: Cambridge University Press. |
[] | Thorne M S, Lay T, Garnero E J, et al. 2007. Seismic imaging of the laterally varying D″region beneath the Cocos Plate[J]. Geophysical Journal International, 170(2): 635–648. DOI:10.1111/gji.2007.170.issue-2 |
[] | Tsuboi S, Geller R J, Morris S P. 1985. Partial derivatives of the eigenfrequencies of a laterally heterogeneous earth model[J]. Geophysical Research Letters, 12(12): 817–820. DOI:10.1029/GL012i012p00817 |
[] | Yan Z Z, Zhang H, Yang C C, et al. 2008. An initial study of the numerical simulation of the Earth's free oscillations process excited by earthquake[J]. Advances in Earth Science (in Chinese), 23(10): 1020–1026. DOI:10.11867/j.issn.1001-8166.2008.10.1020 |
[] | 陈世仲, 李小凡, 汪文帅. 2016. 地球自由振荡的保结构模拟[J]. 地球物理学报, 59(5): 1685–1695. DOI:10.6038/cjg20160513 |
[] | 方俊. 1990a. 地球自由振荡线性反演中的参数核(Ⅰ)[J]. 中国科学B辑, 20(2): 202–207. |
[] | 方俊. 1990b. 地球自由振荡线性反演中的参数核(Ⅱ)[J]. 中国科学B辑, 20(3): 329–336. |
[] | 倪四道. 2007. 超低速区研究进展与展望[J]. 中国科学技术大学学报, 37(8): 820–829. |
[] | 严珍珍, 张怀, 杨长春, 等. 2008. 地震激发地球自由振荡过程的数值模拟初步探索[J]. 地球科学进展, 23(10): 1020–1026. DOI:10.11867/j.issn.1001-8166.2008.10.1020 |