图1 全球大地水准面观测模型EGM2008[8]
Fig. 1 Global geoid observation model EGM2008[8]
北京大学学报(自然科学版) 第61卷 第2期 2025年3月
Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 61, No. 2 (Mar. 2025)
doi: 10.13209/j.0479-8023.2025.011
国家自然科学基金(92155204)资助
收稿日期: 2024–02–22;
修回日期: 2024–04–17
摘要 针对特提斯构造域大尺度重力场响应问题, 通过对全球大地水准面的球谐谱和核函数分析, 将全球地震层析成像模型作为地幔结构, 进行大地水准面正演计算, 探究由特提斯洋动力过程产生的地幔结构与当今特定球谐阶数的特提斯域大地水准面异常之间可能的联系。全球大地水准面观测的球谐谱分析结果表明, 特提斯构造域的大地水准面异常具有球谐 7 阶的波长尺度。通过全球大地水准面核函数的分析, 初步约束特提斯构造域大地水准面异常的密度源深度。最后, 通过全球地球动力学数值模拟计算, 测试并验证地幔 425～1000km 深度处的球谐 7 阶密度异常与特提斯构造域大地水准面异常在印度尼西亚、缅甸、青藏高原及扎格罗斯地区具有很好的对应性, 认为目前观测到的球谐 7 阶大地水准面异常指示位于地幔 425～1000km 深度的特提斯洋中–东段的板片残余信号。
关键词 大地水准面; 特提斯构造域; 球谐谱; 核函数; 地球动力学数值模拟
特提斯构造域从新几内亚岛起始, 经过爪哇和苏门答腊岛, 向北延申到缅甸地区, 进入青藏高原南缘, 再经由伊朗高原, 向西延伸至地中海, 最终绵延到欧洲西部, 总长度达 15000km。这一构造域是当今地球上环太平洋构造域以外规模最大的构造域, 具有全球大陆最丰富的地质现象, 同时蕴藏着丰富的矿产和油气资源[1–3]。特提斯构造域是多期特提斯洋闭合的产物, 其复杂的构造特征也源于多期古大洋的演化和最终闭合。已经有大量基于地质学和地震学观测结果的研究, 讨论不同期次特提斯洋的演化和闭合的动力学过程以及这些过程对现今地幔结构产生的影响[4–7]。然而, 关于特提斯构造域引起的地球重力场变化, 学界目前关注较少。
本文从反映地球重力位起伏的大地水准面观测结果出发, 结合地震层析成像对地幔结构的约束以及瞬时 Stokes 流的地球动力学数值模拟计算, 约束特提斯构造域的重力场响应及其质量异常的来源, 进一步理解特提斯洋演化的动力过程及其对现今地球物理场的影响。
大地水准面指地球重力场中与处于自由静止状态的平均海水面重合或最接近的重力等位面。根据目前常用的大地水准面观测模型 EGM2008 (图 1)[8], 大地水准面的起伏在西伯利亚、印度洋北部以及南极洲附近呈现负值, 其中印度洋北部的负值区域表现得尤其显著, 该区域也称为印度洋大地水准面低值区(Indian Ocean geoid low, IOGL)。这些负异常形成的原因是近年来大地水准面研究的热点问题, 前人对此已有充分的讨论[9–12]。然而, 关于位于青藏高原区域, 将 IOGL 和西伯利亚这两个低值异常区隔开的相对高值区的研究相对欠缺。由于这个相对高值区域正好与当今特提斯构造域在位置上有一些重合, 我们猜测它与特提斯洋的动力演化可能存在联系, 所以本研究聚焦这部分大地水准面正异常产生的原因。
为了分析上述区域大地水准面异常产生的原因,我们首先对该区域信号的球谐成分进行分析。分别在青藏高原区域内采集若干位点(图 2(a)中白色圆点), 通过球谐变换, 计算每一样本点中每个球谐项的分量值, 得到每一个阶次对总信号的贡献强度。对于这些样本点,
图1 全球大地水准面观测模型EGM2008[8]
Fig. 1 Global geoid observation model EGM2008[8]
(1)
其中, θ 和f分别为样本点的纬度和经度, h 为样本点的大地水准面异常幅度, l 为最高展开到 20 阶的球谐阶数, m 为球谐次数, alm 为球谐系数的余弦项, blm为正弦项, plm 为归一化的连带 Legendre 多项式。为了分析每个单一阶次分量对上述区域的大地水准面贡献, 我们对这些样本点的每个单一阶次的大地水准面值进行平均。由于采样点按照 2°×2°的网格选取, 其空间分布密度与纬度相关(与 sinθ 成正比), 因此在计算平均值时, 我们用 sinθ 对每个采样点进行加权。l 阶大地水准面平均值用以下公式表示:
图2 青藏高原高值异常的采样点(a)及球谐成分分析结果(b)
Fig. 2 Sampling points of high anomalies in Tibetan Plateau (a) and results of spherical harmonic analysis (b)
(2)
其中, θi 和 fi 为样本点 i 的地理坐标, N 为样本点数目。由于青藏高原的大地水准面异常是两个主要大地水准面低值区之间的一个相对高值区域, 所以我们主要关注采样点各个阶次对大地水准面的正贡献, 忽略所有负贡献的阶次。我们得到归一化的单一阶次正向大地水准面贡献比值:
(3)
比较每个阶次的 Pl, 可以确定青藏高原大地水准面异常的主要波长。图 2(b)展示每个阶次的 Pl 球谐贡献谱, 可见球谐 7 阶信号对青藏高原区域正值异常成分的贡献最大。为了探究大地水准面的球谐 7 阶信号, 我们对 EGM08 模型进行球谐分解, 仅留下球谐 7 阶的信号(图 3(a))。
图 3(a)中, 球谐 7 阶成分的分布形态正如我们的预期, 恰好对应青藏高原的高值异常。值得关注的是, 我们发现球谐 7 阶信号中存在一个正好对应特提斯构造域的高值异常条带。我们将试图探究球谐 7 阶信号中的高值异常与特提斯动力过程之间可能的联系。为了理解图 3(a)所示球谐 7 阶异常, 下面对它们的来源进行分析。
(b)中水平虚线表示地幔转换带深度范围(410～660km)
图3 EGM08 大地水准面观测模型的球谐 7 阶信号(a)、大地水准面核函数(b)和计算核函数及全球大地水准面所使用的一维黏度剖面(c)
Fig. 3 Spherical harmonic degree 7 component of the EGM08 geoid observation model (a), geoid kernel (b), and one-dimensional viscosity profile used for calculating kernel function and global geoid modelling (c)
我们借助图 3(b)所示大地水准面核函数来分析在全球大地水准面观测中产生球谐 7 阶的密度异常。如果地幔的密度结构已知, 对于一个径向黏度剖面(图 3(c)), 可以使用传播矩阵方法[13–14], 分别计算每个球谐阶数和次数的大地水准面。也可以将每个深度的密度展开系数与深度相关的大地水准面核函数相乘, 然后在深度上进行积分来计算每个球谐阶数和次数的大地水准面[15–16]。大地水准面核函数表示地幔密度异常对地表大地水准面的响应, 包括地幔中的密度异常、由密度引起的地表地形的质量异常以及由密度引起的核–幔边界地形的质量异常 3 个方面的贡献[7]。在特定深度 d 和给定阶数 l 处, 当大地水准面核函数值大于 0 时, 表明深度 d 处具有阶数 l 的正(负)密度异常对地表具有阶数l 的大地水准面有正(负)贡献。相反地, 如果大地水准面核函数值小于 0, 则密度异常对地表大地水准面有相反的贡献。据此, 我们可以对上述球谐 7 阶大地水准面异常区进行如下分析: 球谐 7 阶核函数主要在 425～1500km 区间呈现正值, 青藏高原大地水准面分隔区显现的强球谐 7 阶高值异常信号可能来自这一深度区间的正密度异常(或者冷异常、地震波正异常), 对于较浅的深度区间则相反。
为了探究真实地球内部的相关密度结构, 我们统计 18 个不同的全球地震层析成像模型[17], 在不同的深度统计全球每一点出现地震波波速正异常的频次, 绘制地震波速正异常频次图(图 4(a)～(c)), 其中每一点由冷色到暖色代表在该点支持波速正异常的层析成像模型数由少到多。越偏暖色, 则代表该地区在该深度越有可能呈现波速正异常。从图 4 可以看出, 在 670 和 1000km 深度, 多数地震层析成像模型都认为在印度尼西亚和马来西亚区域存在大量地震波速正异常; 在 1500km 深度, 多数地震层析成像模型支持在印度大陆及青藏高原地区存在大量地震波速正异常。我们从图 4(d)中发现, 这两个区域刚好对应大地水准面观测中球谐 7 阶信号高值区, 与球谐 7 阶的大地水准面正异常信号来自这些深度的密度正异常的推测相符。
我们继续猜测上述地幔内部密度异常来自怎样的–地球动力学过程。关于在地震层析成像结果中显示的印度尼西亚、马来西亚、印度大陆及青藏高原地区的波速正异常, 我们认为它们可能与新特提斯洋闭合产生的俯冲作用相关[18]。如图 3(a)所示, 大地水准面观测中的球谐 7 阶正异常信号恰好与现今特提斯构造域的中段和东段位置重合。可以合理地猜想, 特提斯洋中段和东段闭合产生的洋陆俯冲板片与相应区域地震层析成像中的波速正异常相对应, 而这些区域的地震层析成像波速正异常又有可能是地表产生大地水准面球谐 7 阶正异常信号的原因。针对位于欧洲部分的特提斯西段区域, 由于俯冲距离较短, 且在较早时期已结束俯冲, 完成闭合, 其俯冲物质可能存在于更深的地幔当中, 因此在较高的球谐阶数中难以显现。为了验证上述猜想, 我们需要通过定量计算的方法, 将特提斯构造域的演化历史、地幔密度结构及大地水准面的响应信号联系起来, 确定这三者之间的关系。
大地水准面是地球的一个重力等位面, 大地水准面起伏是大地水准面相对于参考面的差距。大地水准面起伏度j服从泊松方程:
(4)
其中, δρtotal 为对大地水准面有贡献的所有密度异常, 其来源包括地球内部密度不均匀分布(δρbuo)、地表地形起伏度(s)(产生密度差 δρt)和核幔边界地形起伏(b)(产生密度异常 δρb)三部分。将 δρtotal, s和 b进行球谐展开:
(5)
(6)
(7)
其中, r为半径, θ为余纬, f为经度, l为球谐阶数, m为球谐次数, 下角标 sin 和 cos 分别表示球谐阶数中的正弦项和余弦项。对式(4)在球谐域中求解, 可得如下公式[19]:
图4 670, 1000 和 1500km 深度波速快异常的“vote map” ((a)～(c))与 EGM08 模型中球谐 7 阶信号(d)的对比
Fig. 4 Comparison of “vote maps” of fast velocity anomalies at depths of 670, 1000 and 1500 km ((a)–(c)) with the spherical harmonic degree 7 signal of the EGM08 model (d)
(8)
式中, r'是源的半径, r 是需要计算位场处距地心的距离, R为地球半径, RCMB 为核幔边界半径, Δρb 和 Δρt分别为核幔边界密度差及地表界面密度差, r<=min (r', r), r>=max(r', r)。
式(8)中的 s 和 b 服从流体力学的质量守恒和动量守恒方程, 可以通过解算斯托克斯流问题得到。在将地幔视为不可压缩流体的情况下, 有
(9)
(10)
式中, u 表示速度, p表示动压, η 和 ρ0 分别表示地幔的黏度和密度, g0表示重力加速度, er表示径向的单位矢量, δρ和δg分别表示密度和径向重力加速度的扰动。
通常使用离散化数值方法解算斯托克斯流问题的数值解。在使用数值代码解算斯托克斯流问题的时候, 我们考虑了动量方程中自重力项(ρ0δg)的影响。它描述地球内部物质分布不均时, 其中一个点重力加速度的变化。根据前人的经验, 自重力效应会在较大程度上影响计算得到的地表地形、核幔边界地形以及大地水准面, 所以我们在求解式(10)时考虑了自重力的影响。此外, 我们并不知道地球内部具有的密度结构(即 δρ), 它是式(10)中求解斯托克斯流的一个条件。我们将通过全球地震层析成像模型的结果获得地球内部的密度结构。
我们使用地幔动力学模拟程序 CitcomS[20]计算大地水准面异常。CitcomS 是一个基于有限元计算方法的程序, 可以在给定密度和黏度分布的情况下求解三维的标准热对流方程。如前所述, 大地水准面的解算构建的是瞬时流动模型, 即只求解质量和动量的守恒方程, 不存在时间相关性。我们通过将地震层析成像得到的波速异常转换为密度异常来驱动地幔的流动。我们的模型使用 65×65×65 个节点, 平均水平分辨率为 1°×1°, 使用自由滑移边界条件。对于 S 波波速模型, 我们使用 0.25 的波速密度转换系数(dlnρ/dlnV)[21–22]。本文使用的地震层析成像模型是 GYPSUM-S[23], 它是通过体波走时和地球动力学观测结果联合反演得到的。大地水准面对径向黏性结构非常敏感, 我们采用的初始一维黏度剖面如图 3(c)所示, 这是一个通过板块驱动的三维地幔对流演化模型得到的随深度变化的一维黏度剖面[24]。此外, 我们还考虑了由 η=η0⋅exp(E(T0−T))给出的黏度–温度关系[25], 其中 η0 是特定层的参考黏度, T0 和T 是无量纲的参考温度和实际温度, E 为活化能, 控制黏度对温度的依赖强度。热膨胀率设定为 3×10−5 K−1, 用于将相对温度异常转换为密度异常。
我们使用 CitcomS 计算包括自重力影响的球谐20 阶的大地水准面, 然后将计算得到的大地水准面与 EGM08 模型观测到的大地水准面进行比较。为了评估计算结果与观测结果的相关性, 我们定义球谐阶数的相关系数:
(11)
(12)
式中, H 为计算得到的大地水准面球谐分量, h 为大地水准面观测的球谐分量。计算得到的全球大地水准面结果如图 5(a)所示, 其与观测得到的经过校正的大地水准面(图 1)的相关系数(C2～20)达到 0.8034。从形态上来看, 计算结果基本上复现现今大地水准面观测当中主要的信号特征。
(a)利用全球地震层析成像模型计算得到的全球大地水准面, 相关系数 C2～20=0.8034; (b)EGM08 大地水准面观测模型的球谐 7 阶信号; (c)利用 GYPSUM-S 层析成像模型计算得到的大地水准面的球谐 7 阶成分, 相关系数 C7=0.7229; (d)GYPSUM-S 模型425～1500km 深度的球谐 7 阶成分驱动产生的大地水准面, 相关系数 C7=0.6650; (e)GYPSUM-S 模型 500km 深度波速异常的球谐 7阶成分; (f)GYPSUM-S 模型 1000km 深度波速异常的球谐 7 阶成分; (g)滤去 425～1500km 深度球谐 7 阶成分后驱动产生的大地水准面, 相关系数 C2～20=0.8003
图5 大地水准面模拟计算结果
Fig. 5 Results of geoid modelling
为了验证本文第 2 节中提出的猜想, 我们对通过地震层析成像计算得到的全球大地水准面计算结果进行球谐分析, 提取其中球谐 7 阶的信号成分(图5(c))。对比计算得到的大地水准面球谐 7阶成分(图5(b))发现, 用层析成像模型计算的结果在特提斯构造域中–西段的信号较明显, 且在全球范围内的相关系数(C7)达到 0.7229, 说明我们用计算的方法复现了第 2 节中提出的大地水准面球谐 7 阶信号。
为了研究计算得到的球谐 7 阶大地水准面异常在地幔当中的密度源, 我们单独提取 GYPSUM-S 模型在 7 阶核函数为正值的深度区间(425～1500km)的球谐 7 阶成分, 在仅保留这些深度的 7 阶成分的地幔结构下, 计算其大地水准面的响应。针对地幔当中的球谐 7 阶结构, 从图 5(e)和(f)显示的层析成像结果可见, GYPSUM-S 模型在 500～1000km 深度的球谐 7 阶正波速异常结构也与特提斯造山带中–西段的位置具有一定的对应性, 说明地震层析成像揭示的该深度区间可能的特提斯俯冲板片残余具有球谐 7 阶的波长特性。这些球谐 7 阶的波速成分可能在一定程度上地代表地幔中的特提斯信号。图 5(d)显示由这些可能代表特提斯俯冲信号的地幔球谐 7阶结构引起的大地水准面响应, 该计算结果在全球范围的相关系数(C7)达到 0.6650。从图 5(d)中可以发现, 在特提斯构造域中–西段的特征区域(包含印度尼西亚、缅甸、青藏高原及扎格罗斯地区), 计算结果均显示大地水准面的正值贡献。这些区域的正值异常也与 EGM08 模型观测及全地幔计算得到的大地水准面结果一致。此外, 当 425～1500km 深度的球谐 7 阶成分被滤去后, 特提斯造山带中–西段大地水准面的计算结果大幅降低(图 5(g)), 进一步说明地表观测到的位于特提斯构造域中–西段的球谐 7 阶大地水准面高值主要是由地幔 425～1000km深度可能指示特提斯中–西段俯冲残余的球谐 7 阶高密度结构驱动产生的。
本文通过对全球大地水准面观测结果进行球谐分析, 发现全球大地水准面球谐 7 阶成分中的正值异常与特提斯构造域的中–西段具有很强的对应关系, 我们猜测这样的正值异常是由特提斯演化闭合过程中保留在地幔中的球谐 7 阶波长的俯冲物质引起的。通过对 18 个全球地震层析成像模型的统计分析发现, 在 1500km 深度以上的波速结构均显示位于特提斯构造域中–西段区域的波速正异常, 这些异常与地表观测到的球谐 7 阶大地水准面异常在位置和尺度上相对应。为了验证特提斯域地幔中的球谐 7 阶结构是特提斯域球谐 7 阶大地水准面的密度来源, 本文借助全球地球动力学数值模拟计算, 用地震层析成像模型 GYPSUM-S 作为地幔密度模型, 通过解算 Stokes 流问题, 正演计算得到与观测结果有良好对应性的全球大地水准面形态。计算得到的全球大地水准面的球谐 7 阶成分与观测结果一致, 均在现今特提斯构造域中–西段区域显示正值异常。本文进而探究仅由地幔中的球谐 7 阶密度结构驱动得到的大地水准面信号, 发现层析成像中的球谐 7 阶波速信号与特提斯构造域具有对应关系, 可以在一定程度上代表特提斯俯冲残余。本文发现在仅保留可能代表特提斯信号的 425～1500km 深度地幔球谐 7 阶结构的情况下, 驱动得到的大地水准面同样在特提斯构造域的特征区域显示明显的正异常, 说明特提斯俯冲的球谐 7 阶波长的地幔结构是引起特提斯域中西段大地水准面正值的来源。
关于通过地震层析成像获得的特提斯构造域地幔结构的地质解释, 仍然存在讨论的空间。这些观测到的高波速结构是否由特提斯动力过程产生(如青藏高原下方的高密度结构就有可能是由印度板块俯冲或青藏板块拆沉引起), 是由哪一期特提斯洋的打开和闭合引起, 也有较大的不确定性。本文工作是一次初步的探索, 目前还无法将地幔结构与地质过程精准地联系在一起。未来的工作中, 计划使用全球古板块重建的地表速度场模型, 模拟相对真实的特提斯洋演化闭合的动力学过程, 利用动力学计算得到的地幔结构, 动态地追踪特提斯洋不同演化阶段的大地水准面变化。这样, 可以进一步将现今特提斯构造域的地幔结构与地质过程紧密地联系在一起, 直接揭示由特提斯洋的演化产生的大地水准面信号, 进一步验证本文提出的猜想。
参考文献
[1] 吴福元, 万博, 赵亮, 等. 特提斯地球动力学. 岩石学报, 2020, 36(6): 1627–1674
[2] 朱日祥, 赵盼, 赵亮. 新特提斯洋演化与动力过程. 中国科学: 地球科学, 2022, 52(1): 1–25
[3] Sengör A M C. Tectonics of the Tethysides: orogenic collage development in a collisional setting. Annual Review of Earth and Planet Science, 1987, 15: 213–244
[4] Metcalfe I. Gondwana dispersion and Asian accretion: tectonic and palaeogeographic evolution of eastern Te-thys. Journal of Asian Earth Sciences, 2013, 66: 1–33
[5] Becker T W, Faccenna C. Mantle conveyor beneath the Tethyan collisional belt. Earth and Planetary Science Letters, 2011, 310(3/4): 453–461
[6] de Jonge M R, Wortel M J R, Spakman W. Regional scale tectonic evolution and the seismic velocity struc-ture of the lithosphere and upper mantle: the Mediter-ranean region. Journal of Geophysical Research, 1994, 99: 12091–12108
[7] Hafkenscheid E, Wortel M J R, Spakman W. Subduc-tion history of the Tethyan region derived from seismic tomography and tectonic reconstructions. Journal of Geophysical Research, 2006, 111(B8): B08401
[8] Pavlis N K, Holmes S A, Kenyon S C, et al. The de-velopment and evaluation of the Earth Gravitational Model 2008 (EGM2008). Journal of Geophysical Re-search: Solid Earth, 2012, 117: B04406
[9] Rao B P, Kumar M R. Seismic evidence for slab grave-yards atop the core mantle boundary beneath the Indian Ocean Geoid Low. Physics of Earth and Planetary Interiors, 2014, 236: 52–59
[10] Reiss A S, Thomas C, van Driel J, et al. A hot midman-tle anomaly in the area of the Indian Ocean Geoid Low. Geophysical Research Letters, 2017, 44 (13): 6702–6711
[11] Steinberger B, Rathnayake S, Kendall E. The Indian Ocean Geoid Low at a plume-slab overpass. Tectono-physics, 2021, 817: 229037
[12] Pal D, Ghosh A. How the Indian Ocean geoid low was formed. Geophysical Research Letters, 2023, 50: e2022GL102694
[13] Hager B H, O’Connell R J. Kinematic models of large-scale flow in the Earth’s mantle. Journal of Geophy-sical Research: Solid Earth, 1979, 84(B3): 1031–1048
[14] Hager B H, R J O’Connell. A simple global model of plate dynamics and mantle convection. Journal of Geo-physical Research, 1981, 86: 4843–4867
[15] Ricard Y, Fleitout L, Froidevaux C. Geoid heights and lithospheric stresses for a dynamic Earth. Annales Geophysicae, 1984, 2: 267–286
[16] Richards M A, Hager B H. Geoid anomalies in a dy-namic Earth. Journal of Geophysical Research, 1984, 89: 5987–6002
[17] Hosseini K, Matthews K J, Sigloch K, et al. SubMachi-ne: web-based tools for exploring seismic tomography and other models of Earth’s deep interior. Geochemi-stry, Geophysics, Geosystems, 2018, 19(5): 1464–1483
[18] van Hinsbergen D J, Lippert P C, Li S, et al. Recon-structing Greater India: paleogeographic, kinematic, and geodynamic perspectives. Tectonophysics, 2019, 760: 69–94
[19] Zhong S, McNamara A, Tan E, et al. A benchmark study on mantle convection in a 3-D spherical shell using CitcomS. Geochemistry, Geophysics, Geosys-tems, 2008, 9: Q10017
[20] Zhong S, Zuber M T, Moresi L, et al. Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection. Journal of Geophysical Research, 2000, 105(B5): 11063–11082
[21] Ghosh A, Becker T W, Zhong S J. Effect of lateral vis-cosity variations on the geoid. Geophysical Research Letters, 2010, 37: L01301
[22] Ghosh A, Becker T W, Humphreys E. Dynamics of the North American continent. Geophysical Journal Inter-national, 2013, 194: 651–669
[23] Simmons N A, Forte A M, Boschi L, et al. GyPSuM: a joint tomographic model of mantle density and seismic wave speeds. Journal of Geophysical Research, 2010, 115: B12310
[24] Liu X, Zhong S. Constraining mantle viscosity struc-ture for a thermochemical mantle using the geoid ob-servation. Geochemistry, Geophysics, Geosystems, 2016, 17: 895–913
[25] Liu, X, Zhong S. The long-wavelength geoid from three-dimensional spherical models of thermal and thermochemical mantle convection. Journal of Geo-physical Research: Solid Earth, 2015, 120(6): 4572–4596
Preliminary Investigation on the Origin of Geoid Anomalies in the Tethys Tectonic Domain
Abstract Regarding the large-scale gravity field response in the Tethyan tectonic domain, this work aims to connect specific wavelength of geoid anomaly with the mantle structure caused by Tethys evolution. To acheive this, spherical harmonic spectrum analysis and kernel analysis on global geoid are employed, and the global geoid is further calculated by forward modeling with the mantle structure obtained by global seismic tomography. Density structure in our geoid calculation is obtained by seismic tomography. Through spherical harmonic analysis of global geoid observation, the geoid anomalies of the Tethyan tectonic domain exhibit a spherical harmonic degree 7 wavelength scale. Analysis of geoid kernels provides initial constraints on the density source depth of Tethyan geoid anomalies. Finally, through global geodynamic numerical modeling, the degree 7 density anomaly at depths of 425–1000 km in the mantle corresponds well with the Tethyan geoid anomalies in the Indonesian, Burmese, Tibetan Plateau, and Zagros regions. Preliminarily, the observed degree 7 geoid anomalies indicate slab remnants in the central and eastern segments of the Tethys at depths of 425–1000 km in the mantle.
Key words geoid; Tethyan tectonic domain; spherical harmonic spectrum; kernel function; geodynamic numerical modelling