crust 1.0 是一个全球的三维地壳模型,其分辨率是 1 度。更老的版本还有分辨率低一些的 crust 5.1 和 crust 2.0。

官方网址:http://igppweb.ucsd.edu/~gabi/crust1.html

模型基本信息

  • 模型将全球划分为 1 度 x1 度的网格,共计 360x180=64800 个单元;
  • 模型分为 9 层,每一层给出 P 和 S 波速度以及密度。这 9 层分别为
    1. 水层;
    2. 冰层;
    3. 上沉积层;
    4. 中沉积层;
    5. 下沉积层;
    6. 上地壳;
    7. 中地壳;
    8. 下地壳;
    9. 地幔
  • 模型中的海水深度和地形数据来自于 NOAA 的 etopo1 数据;etopo1 是 1 分精度的全球地形起伏模型, crust1.0 从该模型中提取海水深度和地形数据,做平均处理得到 1 度单元内的值;
  • 模型中的冰层厚度来源于 etopo1 的 bedrock 和 ice surface 两个模型的差值;
  • 地幔的参数,来自于 LLNL-G3Dv3 模型;
  • 模型中还包含了每个单元的地壳类型信息;地壳类型信息可能主要由地壳年龄决定;

下载

crust1.0 提供了多个压缩包,分别包含了不同的内容:

  1. crust1.0.tar.gz:主压缩包,包含了全球地壳速度和密度模型;
  2. crust1.0-addon:附加包,包含了全球地壳类型模型;
  3. crsthk.xyz:地壳厚度的 XYZ 文件(不含水厚度);
  4. depthtomoho.xyz:Moho 深度的 XYZ 文件(相对于海平面);
  5. sedthk.xyz:沉积层厚度的 XYZ 文件(单位 km);
  6. sedthk-m.xyz:沉积层厚度 XYZ 文件(单位为 m);

crust1.0

crust1.0 中含说明文档 readme 一份,模型文件四个,用于读取模型的 Fortran 代码三个。

模型

模型文件有四个,从模型文件的名字也可以猜出来,四个文件分别是 Vp、Vs、密度和边界的深度。

  • 全球被划分为 1 度 x1 度的网格单元;
  • 网格点定义在单元的中心,即纬度 5-6 度、经度 150-151 度的单元的值位于网格点 (5.5,150.5) 处,也就是 GMT 中所说的像素配准;
  • 全球共计 360*180 个单元;
  • 数据点按照经度优先的准则保存,即在循环时先循环纬度,再循环经度;
  • 第一个点的坐标为 (179.5W,89.5N),第二个点的坐标是 (178.5,89.5),依次类推;
  • 每个模型文件中共有 64800 行,对应 64800 个单元,每行 9 列,对应每一层的值;
  • crust1.bnds
    • 第一列:水层边界的海拔
    • 第二列:水层边界的海拔
    • 第三列:冰层边界的海拔
    • 第四列:上沉积层边界的海拔
    • 第五列:中沉积层边界的海拔
    • 第六列:下沉积层边界的海拔
    • 第七列:上地壳边界的海拔
    • 第八列:中地壳边界的海拔
    • 第九列:下地壳边界的海拔

程序

编译代码

使用如下命令编译:

gfortran getCN1point.f -o getCN1point

其他同理。

getCN1point

getCN1point 用于获取任意经纬度处的地壳模型,具体用法如下例:

$ ./getCN1point
  .... reading all maps ...
 enter center lat, long of desired tile (q to quit)
50 100
#注意纬度在前,经度在后;如果搞反了,第一个值超过了 90,程序也不会报错而是 “正常地” 显示结果
ilat,ilon,crustal type:   41 281
 topography:    1.6400000
  layers: vp,vs,rho,bottom
   1.50   0.00   1.02   1.64
   3.81   1.94   0.92   1.64
   2.50   1.07   2.11   1.54
   0.00   0.00   0.00   1.54
   0.00   0.00   0.00   1.54
   6.10   3.55   2.74 -19.54
   6.30   3.65   2.78 -38.22
   7.00   3.99   2.95 -46.36
 pn,sn,rho-mantle:    7.99   4.44   3.30
 enter center lat, long of desired tile (q to quit)

getCN1maps

getCN1maps 从 4 个模型文件中提取信息,生成多个 Z 文件。

  1. 生成各层的 Vp、Vs、$\rho$、边界深度,计 4*9=36 个文件,文件名 map-vp[n] 代表第 n 层的 Vp,其他类似;
  2. 生成各层的厚度,计 1*8 个文件,文件名类似 map-th[n],第 n 层的厚度由第 n+1 个边界的深度减去第 n 个边界的深度的结果取负值得到;
  3. 生成沉积层厚度 sedthk,由 3-5 层的厚度相加得到;
  4. 地壳厚度 crsthk:冰层 + 沉积层 + 6-8 层厚度

生成的 46 个文件均为 ASCII 格式,只有 Z 值,没有经纬度坐标。可以通过 GMT 的 xyz2grd 命令转换成 GMT 可识别的 netCDF 格式。

GMT4:

xyz2grd crsthk -Rd -I1/1 -Gout.grd -ZTLA -F -V

GMT5:

gmt xyz2grd crsthk -Rd -I1/1 -Gout.grd -ZTLA -r -V

说明:

  1. 使用 -Rd-R-180/180/-90/90 均可,但不可使用 -Rg
  2. 注意 -ZTLA 选项的含义;
  3. GMT5.1.1 的 xyz2grd 存在 Bug,因而该命令仅在 GMT5.1.2 及其之后版本中可用。

getCN1xyz

getCN1maps 生成类似的文件,只是此时的文件为 xyz 文件,每行三列。文件名以 xyz 开头或结尾。XYZ 文件相对来说更易读,因而推荐使用 getCN1xyz 而不是 getCN1maps

将 xyz 文件转换为 GMT 可识别的网格文件,使用 xyz2grd 。注意与上面命令的区别。

GMT 4:

xyz2grd crsthk.xyz -Rg -I1/1 -Gout.grd -F -V

GMT 5:

gmt xyz2grd crsthk.xyz -Rg -I1/1 -Gout.grd -r -V

说明

程序输出的地壳模型还是很让人困惑的,这里用 getCN1point 获得的某一点的模型,并对输出结果做细致地解释。

对于 (100.5E, 50.5N) 来说:

$ ./getCN1point
  .... reading all maps ...
 enter center lat, long of desired tile (q to quit)
50.5 100.5
ilat,ilon,crustal type:   40 281
 topography:    1.80999994
  layers: vp,vs,rho,bottom
   1.50   0.00   1.02   1.81
   3.81   1.94   0.92   1.81
   2.50   1.07   2.11   1.71
   0.00   0.00   0.00   1.71
   0.00   0.00   0.00   1.71
   6.10   3.55   2.74 -18.93
   6.30   3.65   2.78 -37.22
   7.00   3.99   2.95 -45.19
 pn,sn,rho-mantle:    7.96   4.43   3.28

需要注意,第四列给出的是每一层的下边界的海拔。记住这一点,就可以从输出中提取出很多信息:

  1. 这一点的地形为 1.81 km,注意,这里实际上是一度范围内的平均地形;
  2. 水层的下边界深度是 1.81 km,与地形相同,所以水层厚度为零;
  3. 冰层的下边界深度是 1.81 km,与水层的下边界深度相同,所以冰层厚度为零;
  4. 上沉积层的下边界深度是 1.71 km,所以上沉积层厚度为 0.1 km;
  5. 中沉积层和下沉积层厚度均为 0 km;
  6. 上地壳的下边界深度为 18.93 km,算是地形并减去沉积层,上地壳的厚度是 20.64 km;
  7. 中地壳厚度为 18.29 km,下地壳的厚度为 7.97 km;

几个常见的疑问:

  1. 为何中、下沉积层的速度和密度为零?

    因为此处中、下沉积层的厚度为零,即不存在这两层,不存在的东西当然不用给速度和密度了。

  2. 为何水层和冰层的厚度为零,但是却有速度和密度?

    虽然此处水层和冰层的厚度为零,但是因为水和冰的速度和密度在全球范围内是一个常数, 所以虽然这里没有水和冰,还是可以给一个正确的速度和密度的。不像沉积层,不同地方的 速度和密度差很大。这一点可以通过查看全球水层和冰层的速度和密度极值来验证。

绘图示例

#!/bin/bash
gmt grd2cpt out.grd -Cpolar > out.cpt
gmt grdimage out.grd -Rd -JN6i -Bx60 -By30 -Cout.cpt -V -K > a.ps
gmt pscoast -R -J -W0.1p -O >> a.ps

没有认真选择 cpt 文件,看上去效果不好,从细节上看,数据的转换是没有问题的。

修订历史

  • 2013-10-03:初稿;
  • 2014-06-10:加入了 GMT5 的命令;
  • 2015-05-14:重新整理了文章结构,并对模型做更细致解释;