岩石侧压系数
Sheorey (1994) 基于一维层状模型，考虑地球的曲率、壳幔物质属性、温度场和与热相关的其他性质，提出的一种计算地壳初始应力的方法:
k = 0.25 + 7E(0.001 + 1/z)
下图黑线是Tobias Hergert在2011年一篇论文中使用该公式画出的曲线：
k-h
- k 为岩石的侧压系数，是平均水平应力(SH + Sh)/2 与垂向应力SV 的比值(SH 和Sh 分别为最大和最小水平应力)；
- E为岩石的杨氏模量（GPa），这里采用模型的平均岩石模量E=96 GPa, 该平均值依据CRUST1.0 计算得到；
- z 为距离地表的深度，公式中采用以米为单位
GMT绘图
当E=95Ga，深度z以km为单位时，公式可换算如下：
z = 0.665/(k-0.915)
#!/bin/bash
data=k-values.csv
gmt begin
gmt set FONT_ANNOT_PRIMARY=8p,4
gmt set FONT_LABEL=9p,4
gmt set MAP_TICK_LENGTH -2p
gmt set MAP_FRAME_PEN 0.5p
gmt set MAP_GRID_PEN_PRIMARY 0.3p,gray,-
gmt figure k-values pdf A0.2c,E800
gmt basemap -R0/3/0/8 -JX10c/-8c -Bxag0.5+l"k\075\050S@-H@-\053S@-h@-\051\0572S@-v@-" -Bya2fg1+l"depth(km)" -BWNrb
# 绘制实测应力点
awk -F"," '{if(NR<=25) print $1,$2}' $data | gmt plot -Sc0.15c -W0.5p,blue -lKTB+jBR+gwhite+f8p+p0.5p,gray+o0.5p
awk -F"," '{if(NR>=26 && NR<=30)print $1,$2}' $data | gmt plot -Sd0.2c -W0.5p,red -lSAFOD
awk -F"," '{if(NR>=27) print $1, $2}' $data | gmt plot -Sx0.2c -W0.2p,gray -l"Global Dataset"
# 绘制测压曲线
# k = 0.25 + 7E(0.001 + 1/z) E单位为GPa,Z单位为m;
# z = 0.665/(k-0.915) ; E=95GPa, z单位为km
gmt math -T0.915001/3/0.01 T 0.915 SUB INV 0.665 MUL = | gmt plot -W1p -l"E\07595GPa"
# k=1/3
gmt plot -W0.6p,- -l"k\0751/3" << EOF
0.333333 8
0.333333 0
EOF
gmt end
k-z
