潮汐勒夫数与荷载勒夫数的计算


潮汐勒夫数(Tidal Love Numbers, TLNs)反映了地球对潮汐的形变响应,通常用 (h, l, k) 来表示起因(潮汐)与响应(垂直位移、水平位移、引力势扰动)的比例系数。与之相似,荷载勒夫数(Load Love Numbers, LLNs)或荷载形变系数(Load Deformation Coefficients, LDCs)反映了地球对表面荷载的形变响应,通常用 ($h’$, $l’$, $k’$) 来表示起因(表面荷载)与响应(垂直位移、水平位移、引力势扰动)的比例系数。勒夫数与地球的分层结构(壳幔核的分层数目与厚度)以及每层的物质特性(密度、刚性、粘滞性)密切相关,不同的地球模型对应不同的勒夫数。对于球对称弹性不可压缩的自引力(Spherically symmetric, Elastic, Incompressible, and Self-Gravitating, SEISG)地球模型而言,弹性勒夫数(Elastic Love Numbers)只与地球的密度和剪切模量(刚性)有关;对于球对称粘弹性不可压缩的自引力(Spherically symmetric, Viscoelastic, Incompressible, and Self-Gravitating, SVISG)地球模型而言,粘弹性勒夫数(Viscoelastic Love Numbers)不仅与地球的密度和剪切模量(刚性)有关,还与地球的粘滞性(Viscosity)有关。

PREM 地球模型

PREM(Preliminary Reference Earth Model) 地球模型以分段多项式的形式定义了密度、S 波波速以及 P 波波速随地球半径的变化,通过公式 $$\left \{ \begin{array} { l } { \lambda = \rho \left( V _ { P } ^ { 2 } – 2 V _ { S } ^ { 2 } \right) } \\ { \mu = \rho V _ { S } ^ { 2 } } \\ { g ( r ) = \frac { 4 \pi G } { r ^ { 2 } } \int _ { 0 } ^ { r } \rho \left( r ^ { \prime } \right) {r ^ \prime}^2 dr ^ { \prime } } \\ { p ( r ) = \int _ { r } ^ { a } \rho \left( r ^ { \prime } \right) g \left( r ^ { \prime } \right) dr ^ { \prime } } \end{array} \right.$$ 可计算第一拉莫数(Lamé parameters) $\lambda$ 和第二拉莫数 $\mu$(在流体动力学中称为动力学粘滞性,在弹性力学中称为剪切模量)、重力加速度以及压力随地球半径的变化。根据 PREM 模型的分段多项式,下面的 Python 代码定义了函数 $\rho(r)$、$V_S(r)$ 以及 $V_P(r)$

计算密度、S 波波速以及 P 波波速在不同深度的值并作图

输出的图像为

上图显示密度、S 波波速、P 波波速在核幔边界(r=3480km)发生了明显跃变。

用 ALMA 计算潮汐勒夫数与荷载勒夫数

潮汐勒夫数与荷载勒夫数的理论计算主要涉及到质量守恒方程(连续性方程)、动量守恒方程、角动量守恒方程(Liouville 方程)、本构方程(应力与应变的关系)、拉普拉斯方程(引力势的球谐展开)的解算以及边界条件的设定,详细的数值算法可参见 Practical Numerical Computation of Love Numbers & Applications。ALMA 是一套专门用来计算潮汐勒夫数与荷载勒夫数的 Fortran 程序(ALMA 中表 1),可联系 Giorgio Spada(Email: giorgio.spada@gmail.com) 获得,或点击这里进行下载。与 Spada 的源程序有所不同,这里对其做了以下扩展:

  • 在 visco.f90 中添加了 ‘kv=5’ 模式使其能够计算 PREM 地球模型的弹性勒夫数;
  • 修改了 love_numbers.f90 程序,使其在 ‘imode=3’ 的模式下可计算服从线性流变(Linear rheology)的多层地幔地球模型的勒夫数,而不仅仅是计算单一层均匀地幔地球模型的勒夫数;

使用 ALMA 计算潮汐勒夫数与荷载勒夫数,首先要配置文件 alma.inc,以计算 PREM 地球模型的弹性荷载勒夫数为例,设置地球模型为 PREM: imode=1、岩石圈厚度为 120km: lth=120.e3,弹性地幔: kv=5,大陆性地壳: ioc=0,地幔分层数为 11: nal=11,荷载响应: iload=1,荷载历史为 Heaviside: th='h',勒夫数的阶数为 $1\sim 5$、步长为 1: l1=1, l2=5, ls=1,有效数字位数为 64 位(若勒夫数的阶数大于 300 阶,则有效数字位数需设置为 128 位): nsd=64,时间点的采样数为 11: p=10,时间范围为 $10^{-3}\sim 10^3$kyrs: m1=-3., m2=+3.。ALMA 中地球的核(r ≤ 3480 km)默认为无粘滞性(Inviscid),计算 PREM 地球模型分层参数的策略为

  • 将地幔分为等厚度的多层;
  • 用体积平均(volume-averaging)法计算每一层的密度、剪切模量、粘滞性等参数,即 $\overline{\mu} = \frac{\int ^b_a \mu(r)r^2dr}{\int ^b_a r^2dr}$。

配置文件 alma.inc 的样板如下

进入 ALMA 主目录,终端输入 bash FMLIB_compile.sh,屏幕上打印出以下信息

终端输入 bash alma.sh,屏幕上打印出以下信息

并且主目录会生成多个文件(见 ALMA 表 2,表中的 ‘-‘ 表示串联, ‘|’ 表示并联),其中 alma.log 的样板如下

荷载勒夫数文件:h.dat,l.dat,k.dat,其中 h.dat 的样板如下

l.dat 的样板如下

k.dat 的样板如下

输出文件显示弹性勒夫数不随时间发生变化,这是因为弹性勒夫数的表达式为 $$
\left \{ \begin{array} { l } { h _ { l } } \\ { l _ { l } } \\ { k _ { l } } \end{array} \right \} ( t ) = \left \{ \begin{array} { l } { h _ { l } } \\ { l _ { l } } \\ { k _ { l } } \end{array} \right \} ^ { E } H ( t ) \;, $$ 其中 $l$ 为阶数,$H ( t )$ 为 Heaviside 函数。

以计算 PREM 地球模型的粘弹性荷载勒夫数为例,设置地球模型为 PREM: imode=1、岩石圈厚度为 120km: lth=120.e3,粘弹性地幔(Continuous viscosity profile): kv=1,地幔分层数为 3: nal=3,荷载响应: iload=1,荷载历史为 Heaviside: th='h',勒夫数的阶数为$1\sim 5$、步长为 1: l1=1, l2=5, ls=1,有效数字位数为 64 位(若勒夫数的阶数大于 300 阶,则效数字位数需设置为 128 位): nsd=64,时间点的采样数为 101: p=100,时间范围为 $10^{-3}\sim 10^3$kyrs:m1=-3., m2=+3.,地球的核(r ≤ 3480 km)默认为无粘滞性(Inviscid),地幔每一层的密度、剪切模量、粘滞性等参数实际上为体积平均的 PREM 参数。配置文件 alma.inc 的样板如下

终端输入 bash alma.sh,屏幕上打印出以下信息

对输出的 h.dat、l.dat、k.dat 作图

输出的图像为

图中显示粘弹性勒夫数随时间的演化发生变化,这是因为粘弹性勒夫数的表达式为 $$\left\{ \begin{array} { c } { { h } _ { l }} \\ { { l } _ { l }} \\ {{ k } _ { l } } \end{array} \right\} ( t ) = H(t)\left ( \left\{ \begin{array} { c} { h _ { l } } \\ { l _ { l } } \\ { k _ { l } } \end{array} \right\} ^ { F}+\sum _ { i = 1 } ^ { M } \frac{1}{s_{li}}e ^ { s _ { l i } t } \left\{ \begin{array} { l } { h _ { l i } } \\ { l _ { l i } } \\ { k _ { l i } } \end{array} \right\} ^ { V } \right )\;,$$ 其中 $$\left\{ \begin{array} { c} { h _ { l } } \\ { l _ { l } } \\ { k _ { l } } \end{array} \right\} ^ { F} = \left\{ \begin{array} { c} { h _ { l } } \\ { l _ { l } } \\ {k _ { l } } \end{array} \right\} ^ { E} – \sum _ { i = 1 } ^ { M } \frac{1}{s_{li}} \left\{ \begin{array} { l } { h _ { l i } } \\ { l _ { l i } } \\ { k _ { l i } } \end{array} \right\} ^ { V } $$ 为流体荷载形变系数(Fluid LDCs),M 为与地球的分层结构相关的模数。

若采用用户自己提供的地球模型,那么需要将 imode=1 修改为 imode=2imode=3 。ALMA 默认使用 Maxwell(弹簧与阻尼串联的形式)流变来近似地幔的流变特性,若要使用其他的流变特性,需要在 imode=3 的模式下指定流变的类型(见 ALMA 的表 7 和 表 8): ire=i,其中 $i$ 为 1~9 的整数。以 BRL06 为参考地球模型,计算地幔流变特性为 Burgers 的潮汐勒夫数,设置地球模型为 BRL06: imode=3mod3='./rheological_profiles/BRL06.dat',地幔流变特性: re=7,地幔分层数为 5(需要与 BRL06 模型中地幔分层数保持一致): nal=5,潮汐响应: iload=0,荷载历史为 Heaviside: th='h',勒夫数的阶数为 $2\sim 500$、步长为 5: l1=2, l2=500,ls=5,有效数字位数为 128 位(若勒夫数的阶数大于 300 阶,则效数字位数需设置为 128 位): nsd=128,时间点的采样数为 2: p=3,时间范围为 $10^{-3}\sim 10^3$kyrs: m1=-3., m2=+3.

配置文件 alma.inc 的样板如下

终端输入 bash alma.sh,屏幕上打印出以下信息

对输出的 h.dat、l.dat、k.dat 作图

输出的图像如下

图中的潮汐默认服从 Heaviside 函数,而真实的潮汐并不符合 Heaviside 函数,而是永久潮和一些周期项的叠加。


参考

  • Spada G . ALMA, a Fortran program for computing the viscoelastic Love numbers of a spherically symmetric planet[M]. Pergamon Press, Inc. 2008.

关于 “潮汐勒夫数与荷载勒夫数的计算” 的 2 个意见

  1. Pretty nice post. I just stumbled upon your blog and wished to say that I have truly enjoyed surfing around your blog posts. After all I’ll be subscribing to your rss feed and I hope you write again soon!

  2. Every weekend i used to pay a quick visit this website,
    for the reason that i wish for enjoyment, since this site conations really fastidious funny data too.

发表评论

电子邮件地址不会被公开。 必填项已用*标注