全球地壳模型 crust 1.0

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 提供了多个压缩包，分别包含了不同的内容：

crust1.0.tar.gz：主压缩包，包含了全球地壳速度和密度模型；
crust1.0-addon：附加包，包含了全球地壳类型模型；
crsthk.xyz：地壳厚度的 XYZ 文件（不含水厚度）；
depthtomoho.xyz：Moho 深度的 XYZ 文件（相对于海平面）；
sedthk.xyz：沉积层厚度的 XYZ 文件（单位 km）；
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 文件。

生成各层的 Vp、Vs、$\rho$、边界深度，计 4*9=36 个文件，文件名 map-vp[n] 代表第 n 层的 Vp，其他类似；
生成各层的厚度，计 1*8 个文件，文件名类似 map-th[n]，第 n 层的厚度由第 n+1 个边界的深度减去第 n 个边界的深度的结果取负值得到；
生成沉积层厚度 sedthk，由 3-5 层的厚度相加得到；
地壳厚度 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

说明：

使用 -Rd 或 -R-180/180/-90/90 均可，但不可使用 -Rg ；
注意 -ZTLA 选项的含义；
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.81 km，注意，这里实际上是一度范围内的平均地形；
水层的下边界深度是 1.81 km，与地形相同，所以水层厚度为零；
冰层的下边界深度是 1.81 km，与水层的下边界深度相同，所以冰层厚度为零；
上沉积层的下边界深度是 1.71 km，所以上沉积层厚度为 0.1 km；
中沉积层和下沉积层厚度均为 0 km；
上地壳的下边界深度为 18.93 km，算是地形并减去沉积层，上地壳的厚度是 20.64 km；
中地壳厚度为 18.29 km，下地壳的厚度为 7.97 km；

几个常见的疑问：

为何中、下沉积层的速度和密度为零？
因为此处中、下沉积层的厚度为零，即不存在这两层，不存在的东西当然不用给速度和密度了。
为何水层和冰层的厚度为零，但是却有速度和密度？
虽然此处水层和冰层的厚度为零，但是因为水和冰的速度和密度在全球范围内是一个常数，所以虽然这里没有水和冰，还是可以给一个正确的速度和密度的。不像沉积层，不同地方的速度和密度差很大。这一点可以通过查看全球水层和冰层的速度和密度极值来验证。

绘图示例

#!/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：重新整理了文章结构，并对模型做更细致解释；

文章目录