具有任意分层外核的自引力地球中三维波传播的谱元模拟

总结

本文研究全球尺度上地震波传播的谱元模拟。特别强调了与低频研究相关的两个方面。首先，为了充分考虑自引力的影响，将该方法推广到Cowling近似之外。特别是，地球外部重力场的扰动是通过在球面谐波的基础上投影光谱元素解来处理的。其次，我们提出了流体内部的一个新公式，它允许我们考虑任意密度分层。它基于将位移分解为两个标量势，并导致完全显式的流固耦合策略。详细介绍了该方法的实现，并通过一系列基准测试验证了其准确性。

1引言

它是最近几位作者建立的(查尔朱布2000；Komatitsch&Tromp 2002a公司,b条；卡普德维尔2003；Chaljub 2003年)谱元法（SEM）为在地球三维（3-D）模型中计算合成地震图的问题提供了有效的解决方案。当前大多数谱元研究的目的是将计算推向高频，而在全球尺度上传统使用的方法达到了极限，本文重点关注一些对地震频带较低部分至关重要的物理效应：（i）自重的充分处理和（ii）考虑地球流体区域中任何密度分层的能力。

本文的第一个新颖之处在于纳入了自引力，其影响对于周期大于100 s的地震和重力观测非常重要。所有之前基于SEM的研究都解释了Cowling近似下的重力影响(Cowling 1941年)即忽略地震波对重力场的扰动。做出这种假设的主要原因在于问题的内在困难。考虑到自引力的全部影响，需要求解在整个空间上定义的扰动引力势的泊松方程。与球面谐波法不同，使用基于网格的方法（如SEM）并不能为解决外部问题提供自然框架。无界域中基于网格的近似首先通过限制计算域进行，然后通过在截断边界上施加适当的条件进行。根据人工边界条件（ABC）是否是局部的，会出现不同的方法。基于局部ABC的方法具有计算成本低、对任意几何体有效的优点。例如无限元法（例如。贝特斯1992；Gerdes&Demkowicz 1996年)其中，外部解决方案的行为在径向上强制执行。基于非本地ABC的第二类方法并不通用，因为它们通常需要了解外部问题的分析性或半分析性解决方案。它们在精度方面通常具有非常吸引人的特性，但仅限于简单的（通常是球形的）几何图形。非局部ABC可以在Dirichlet to‐Neumann的严格框架内实施到有限元方法中(日期N)操作员（例如。Givoli 1992年). 这是我们在这里保留的方法。这个日期N适合我们问题的算子依赖于地球外拉普拉斯方程解的球面调和分解。与卡普德维尔（2003）为了将与时间相关的谱元计算与频域中的模态解耦合，我们的日期该算子的推导要简单得多，因为它适用于静态问题。泊松-拉普拉斯方程的谱元离散化产生了一个对称的代数系统，该系统必须在每个时间步长进行反演，以获得引力势的扰动。实际上，这是通过迭代共轭梯度法来实现的，其预处理对于执行常规计算至关重要。

我们详细考虑的另一个方面是地核流体部分的处理。关于核心动力学，一个特别重要的参数是平方Brunt–VäisäläfrequencyN个²这是流体对密度扰动的局部响应的特征。就一阶而言，核心可以被视为中性分层，即。N个²=0，因为强烈对流区域的大部分预计会出现中性浮力。然而，有地震学证据表明N个²在核心的顶部N个²批量生产（Masters 1979）。为了通用性，我们对堆芯结构的描述不假设浮力频率的剖面。为此，我们引入了流体中波动方程的双势公式，该公式推广了中性浮力公式Komatitsch&Tromp（2002b）和查尔朱布（2003）与这些考虑流体中速度势的研究相反，我们将分解应用于位移场，以获得扰动重力势的自然固体-流体边界条件。与上述研究相反，这种选择的一个吸引人的结果是产生了一种完全明确的固-液耦合策略。最后请注意，我们的公式接近于Wu&Rochester（1990）在核心动力学研究中，相对于流体区域中的未知数而言，这是最优的。

本文的其余部分组织如下：第2节，我们回顾了自引力地球中强形式和弱形式的运动方程。我们介绍了流体区域位移场的双势分解，并定义了日期N允许我们在有限域内处理方程的运算符。在第3节，我们回顾了空间谱元近似的原理，并详细介绍了我们的显式时间推进算法。最后，数值结果如所示第4节用于验证方法实现的一组球对称模型。

自引力地球中的2波方程

在本节中，我们回顾了波动方程的强形式和弱形式，它是通过围绕非旋转、静水压预应力平衡状态的一阶拉格朗日摄动获得的。参考读者Komatitsch和Tromp（2002年b）用于旋转壳体的延伸。在整篇论文中，地球用？表示，其外边界用？表示_⊕.的固体（或流体）部分称为_S公司（分别）_F类)，所有固液界面的集合表示为∑_旧金山.无论何时考虑地形或椭圆度将表示半径为的球b条包含非球面地球（即。)和将代表其球形边界.

2.1坚固型

2.2边界条件

位移、牵引力和MRP的整套边界条件见Dahlen&Tromp公司（1998年，第104页）。这里我们回顾一下这些涉及MRP或涉及固液界面的边界条件。

2.3弱模板

2.4 日期N操作人员

3数值近似

本节讨论自引力地球中波动方程的数值近似，我们分两步实现。首先，将SEM应用于空间域中方程的弱形式。然后使用有限差分格式及时推进系统。

为了简洁起见，尽可能避免方法的细节。参考读者Komatitsch&Vilotte（1998）和至Komatitsch&Tromp（1999）用于弹性波动方程的SEM的一般描述，以及查尔朱布（2000）,Komatitsch&Tromp（2002a,b）和查尔朱布（2003）将其扩展到全球地震学，包括在具有分布式内存的现代计算机上并行实现。

3.1空间离散化

3.1.1六面体网格

第一个离散化步骤包括将球形地球分解为一组不重叠的六面体元素。该过程详见查尔朱布（2003），其中引入了不符合要求的接口，以避免随着深度对网格进行过度细化。这种策略允许在空间上对网格进行细化（或粗化），复杂性与跨界面不匹配的元素之间的连续性要求有关。为了简单起见，本文仅限于球面几何一致网格的情况，如图6注意，考虑到地球的椭圆形状或考虑到表面地形，在自重情况下需要将网格延伸到人工边界之外.

3.1.2光谱元素法

基于球体的三维铺设，MRP（ψ）以及固体中的位移(u个)用连续张量多项式逼近流体中的势（χ和ξ）。注意，流体区域内法向位移的连续性在弱形式下自然满足(19)和(20).

每个谱元上使用的基多项式被定义为配置点的形状函数。SEM的特殊性之一是，配置点是所谓的高斯-洛巴托-勒让德点，即用于评估弱形式方程中积分的相同点。这种选择的一个结果是L（左）²标量积是对角的，这是一个允许我们设计显式时间方案的属性（参见例如。Komatitsch&Vilotte 1998年；Komatitsch&Tromp 1999年).

3.2时间演变

让我们强调，流体和固体区域之间的耦合不需要迭代等式（31）和(32)如果使用速度势公式（例如。科马提奇2000；查尔朱布2003). 这种吸引人的特性源于潜在的分解(8)应用于位移，这是Newmark方案中的显式变量。

当考虑到自引力的全部影响时，前面的评论仍然有效。从位移场计算MRP确实是明确的，因为它不涉及任何时间导数Ẋ或不用说，这项任务成本很高，因为它需要对对称的病态矩阵进行形式化反演P（P）（例如。德维尔2002). 实际上，我们解决了等式（29）对于采用共轭梯度（CG）方法的MRP，当残差减少一个待选择的系数ε时，迭代停止。本文没有讨论为Poisson–Laplace解算器构建有效的预条件器的问题，但为了避免性能瓶颈，这无疑是至关重要的。

4数值结果

在本节中，我们通过几个例子证明了我们方法的有效性，可以为这些例子导出参考、半分析的解决方案。首先，对于具有恒定Brunt–Väisälä频率的模型，在Cowling近似下测试双势公式，即不计算MRP。其次，PREM模型中包含了质量再分配的影响(Dziewonski&Anderson 1981年).

4.1双电位制剂的验证

为了定义一些基准来测试我们的公式，我们考虑了PREM模型的无海洋版本，如图1该参考模型进一步受到约束，以适应流体外核中布伦特-Väisälä频率平方的给定剖面。为了继续，我们只需改变P（P）速度，c（c），英寸方程（6）保持密度、梯度和重力加速度不变。请注意，调整Brunt–Väisälä频率的现实方法是修改密度，而不是P（P）速度（参见示例。Wu&Rochester 1993年)因为后者在地球上受到更好的约束。然而，根据P（P）速度剖面简单明了，对于数值验证而言仍然完全可以接受。

图2显示了按照上述步骤构建的三个模型。“N”标签指中性分层的外核（即N个²=0），而标有“S”和“U”的模型分别对应稳定层结和不稳定层结。为了简单起见，我们在模型“S”和“U”中选择了平方Brunt–Väisäläfrequency的值为常数，分别等于N个²= 10⁻⁷拉德²秒⁻²和N个²=−5 × 10⁻⁸拉德²秒⁻²这些值对应于地球无地震振荡反演的极值。请注意N个²PREM的相似性表明，PREM内的震级大约小一个数量级P（P）‐中性分层剖面的速度。

这三个模型都是由浅爆炸点源激发的，其时间依赖性是一个主频的Ricker小波（即高斯函数的二阶导数）（f）₀=1毫赫兹。震源位于赤道深处d日_秒≃61 km，接收器沿赤道放置。图3显示了三个模型中在90°震中距离处记录的纵向位移。使用每个模型的特征模式总和，在Cowling近似下计算轨迹。波形差异说明了地震波对液芯分层的敏感性，并表明模型“U”和“S”构成了双势公式的合适基准。在图4和5将这两个模型得到的谱元结果与几个震中距离的模态解进行了比较。这两种解决方案非常接近，在所考虑的时间间隔内，最大的相对差异小到每密耳几个。

用于进行计算的光谱元素网格如所示图6它由1024个元素组成，其中多项式次数在径向上从2（地壳中）到10（外核中）不等，在切向上保持不变，等于8。网格点的总数为443232，对应于每个波长的点数远大于5，这是获得精确解的经验最小比率（例如。Komatitsch&Vilotte 1998年). 这解释了光谱元素计算和参考解之间的匹配。

4.2整个配方的验证

作为最后一个例子，我们考虑了地球模型的弹性重力响应的计算图1该测试包含了本文中提到的所有困难：液芯的分层是任意的，物理描述包括了自重的全部影响。

模拟的参数与上述略有不同，因为震源主频设置为较低的值（f）₀=0.5 mHz，以最大化质量再分配的影响。因此，光谱元素网格被适配，并且在每个方向上与图6.

为了检查测试对实现自引力的要求是否足够高，我们比较了图7记录的表面纵向位移，包括或不包括重力势的扰动。这两条记录道都是通过正常模式总和计算出来的，并在震中距离为90°的地方记录了大约10小时。相位和振幅的差异表明，Cowling近似在实验的频率范围内是无效的。

最后，将SEM获得的结果与参考溶液进行了比较图8考虑了两种情况，这两种情况对应于关于离散泊松-拉普拉斯CG分辨率的光谱元素解的不同精度等式（29）在第一种情况下，一旦残差减少了三个数量级，CG迭代就会停止，这意味着ε=10⁻³得到的光谱元素解显然不够精确，并且包含一个长期项，似乎打破了离散水平上的能量守恒。为了纠正这种行为，我们考虑第二个测试，其中停止标准设置为ε=10⁻⁵在这种情况下，计算在所考虑的时间间隔上是稳定的，光谱元素解的精度是可以接受的，其与参考解的相对差异小于几/mil。

在上述每种情况下等式（26）已设置为我_最大值=20，基于先验的了解PREM中的色散关系。在任何实际应用中我_最大值将始终小于100，因为在小于100 s的周期内，质量重分布的影响可以忽略不计。低估截断阶数的影响是将特征振荡添加到谱元解中（本文中未显示）。值得注意的是，数值误差的两个可能来源（ε太大或我_最大值太小）导致不同的签名。这提供了两种不同的诊断方法，使我们能够构建具有任意精度的光谱元素解决方案。

5结论

我们已经展示了如何调整SEM来解释与全球地震学相关的两个影响：对自重的充分处理和考虑流体外核中任何密度分层的能力。通过在球对称模型中进行的一系列数值试验，证明了该方法的准确性。结合上述两种效应，我们认为基于SEM的数值模拟将为地球三维模型的弹性重力响应提供新的估计。

致谢

欧共体承认与美国普林斯顿大学地震工作队成员进行了多次讨论，这项工作就是在那里开始的。这份手稿得益于卢多维奇·马格林和亚历山大·福尼尔的仔细阅读。计算是在法国格勒诺布尔天文台的计算密集服务委员会（SCCI）和法国蒙彼利埃国家环境信息中心（CINES）进行的。

工具书类