孟秋, 胡才博, 石耀霖
中国科学院大学地球与行星科学学院 中国科学院计算地球动力学重点实验室, 北京 100049
2019年12月20日 收稿; 2020年2月17日 收修改稿
基金项目: 国家重点研发计划项目（2016YFB0200801）、国家自然科学基金重大项目（41590865）和联合基金（U1839207）资助
通信作者: 胡才博, E-mail: hucb@ucas.ac.cn

摘要: 随着全球气候不断变化，祁连山地区广泛存在的现代冰川，可能会进一步消融乃至完全消失，从而引起祁连山地区的地表变形、地应力及断层地震危险性的时空变化。利用自主开发的平面应变黏弹性有限元程序，计算在冰川不断消融引起的地表载荷和两侧地块的构造挤压作用下，1956-2106年共150 a期间，祁连山地区自地表变形、地应力的时空演化。重点关注祁连山地区典型地质剖面的主要断层的正应力、剪应力及库仑应力随时空的变化，定量评估各个断层地震危险性的变化。结果表明，在祁连山冰川不断消融的情况下，祁连山地区地表抬升速率呈现先快速增加后趋于稳定的冰后回弹现象；研究区域的水平应力、垂直应力及剪应力在冰川载荷卸载和边界构造加载的双重作用发生显著变化；大柴旦-宗务隆山断裂、哈拉湖南山断裂和哈拉湖盆地北缘断裂库仑应力变化持续增大，未来地震危险性增强；阴凹槽断裂和祁连山北缘断裂库仑应力增加很小，反映了在冰川融化影响不大、仅仅在构造应力增加为主的情况下，地震危险性随时间增加是缓慢的。
关键词: 有限元黏弹性地震危险性祁连山冰川
Effects of glacier ablation on surface deformation, stress field, and seismic risk of main faults in Qilian Mountain area
MENG Qiu, HU Caibo, SHI Yaolin
Key Laboratory of Computational Geodynamics of Chinese Academy of Sciences, College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing 100049, China


Abstract: With the continuous development of global climate change, the widespread modern glaciers may further melt or even completely disappear, resulting in the temporal and spatial changes of surface deformation, stress, and fault seismic risk in Qilian Mountains. In this paper, a plane-strain viscoelastic finite element program was developed to calculate the spatial and temporal evolution of surface deformation and ground stress in Qilian Mountains during the 150-year period from 1956 to 2106 under the surface load caused by glacier melting and the continuous tectonic compression. We calculated the normal stress, shear stress, and the change of Coulomb stress on the main faults in Qilian Mountains. The changes of seismic risk on each fault were quantitatively estimated. The surface uplift rate of Qilian Mountains is controlled by the ice rebound. The surface uplift rapidly increase at the beginning and then tends to be stable under the continuous melting of Qilian Mountains glacier. The horizontal stress, vertical stress, and shear stress change obviously under the effects of glacier melting and boundary tectonic loading. The Coulomb stress changes ΔCFS of the DZF, HSF, and HNF continuously increase, and the seismic risk increases in the future. Coulomb stress changes on the YAF and QNF faults are very small, and the seismicity seismic risk changes little. The numerical simulation provides a theoretical reference for understanding the spatial and temporal variation of future earthquake risk in Qilian Mountains under the glacier melting and tectonic loading.
Keywords: finite element modelvisco-elasticityseismic riskQilian Mountainsglacier
冰冻圈是指地球上存在的所有冰体所构成的空间，包括包含冰的云层、地下冰、积雪、水冰、冰川、冰盖等瞬时、短期、季节和永久存在的冰体，最早在1923年由Dobrowolsky提出^[1]。人类对冰川和冰盖在调节气候和地球地貌演化等长期效应和泥石流、洪水等短期效应方面都做了深入的研究^[2]。自工业革命以来，与人类活动密切相关的全球气候变化，正在加速冰川消融。研究现代冰川消融乃至消失对地表变形、地应力和主要断层地震危险性的时空变化，具有理论和现实意义。
中国有48 571个冰川，覆盖了5.18万km²的地区，占世界除南极和格陵兰岛之外冰川总和的7.1%^[3]。青藏高原地区拥有世界上除两极之外数量最多的冰川，很多冰川是内陆河流的源头^[4]，对于水循环和物质平衡具有重要的意义^[5]。本文主要关注祁连山地区(36°30′~39°30′N, 93°30′~103°00′E)的现代冰川。祁连山位于中国西北内陆干旱区，隶属青藏高原东北缘，气候属于温带半干旱区，受大陆性气候和青藏高原气候的影响^[6-8]。祁连山地区由许多NW-SE走向的平行山脉和山谷组成，长度约为800 km，宽度约为300 km，东起乌鞘岭，西止于当金山口，向东南延伸到河西走廊，地形从东北向西南逐渐隆升，最高点为团结峰(5 826 m)^[3]。
祁连山地块是青藏高原向北扩张的前缘部位，受到青藏高原与北部阿拉善块体的挤压，发生强烈变形^[9]，GPS整体呈NE-SW向压缩，祁连山北缘断裂滑动主要运动特征为逆冲兼具右旋走滑，呈现出空间不均匀变形的现象^[10]。祁连山主要断裂带包括祁连山断裂带和阿尔金断裂带，其中祁连山断裂带由多条活动断裂带组成，总体为北西西走向伸展，地震活动频繁^[11]。该区域主要受到挤压作用，同时也会受边界走滑作用的影响，形成左旋逆冲兼走滑的构造特征^[12]。祁连山地震带被深大走滑活动断裂包围，区内发生过多次大地震，包括1920年海原8.5级地震、1927年古浪8级地震等^[13]。
祁连山冰川属于阿尔卑斯冰川，是冰冻圈的重要组成部分，同时也是内陆地区重要的水源，被称为“阿尔卑斯固体水库”，其融水对于中国西北部脆弱的山区生态系统起着至关重要的作用，是中国西北部干旱区的重要水源之一。同时，也加强了中国地区与“丝绸之路”其他国家地区的往来，为中国内陆干旱地区经济的发展、文化的繁荣提供了保障^[14]。祁连山冰川中，小于1 km²覆盖面积的冰川数量最多，1~5 km²的冰川覆盖面积最大，其中58%面积的冰川位于海拔4 800~5 200 m之间的区域^[3]。
随着全球气温的升高，祁连山地区的冰川逐渐消融，雪线不断向后退缩，在1956—2000年间可以达到419 m/a。1956—2010年期间，祁连山冰川总面积减少420.81 km²，约占总面积的20.88%；体积减少21.63 km³，占总体积的20.26%。其中小型冰川的减少是主要原因，目前海拔4 000 m以下冰川已经完全消失，4 500 m以下的冰川跟50 a前相比，消失了一半以上。4 350~5 100 m之间的冰川减小面积占总面积的84.21%。整体看来，冰川消融呈现出东部冰川减少最快、西北部减少较慢的趋势^[3]。
祁连山地区如此巨量的冰川载荷变化可能会引起地壳变形、地应力的时空变化，甚至可能诱发新的地震活动。历史上有很多人类活动或者自然环境变化引起地壳变形、影响地震活动性的先例。Liu等^[15]指出，飓风可能引发活动断裂上的慢地震的发生；Bettinelli指出喜马拉雅GPS速度的垂向分量具有季节性，与恒河平原地表水载荷有关^[16]，并且在冬季小震活动更加频繁^[17]。前人对地下水开采、水库蓄放水、地热开发高压注水、油气注水开采等也进行了大量的研究工作。庞亚瑾等^[18]模拟计算华北地区地下水的开采对地壳应力的影响，结果表明地下水的开采有减缓大地震孕育的作用；程惠红等探讨紫坪铺水库对汶川地震的触发作用，结果表明，水库造成的影响恰好处于是否触发的边缘^[19-20]，周围小震应该与紫坪铺蓄水有直接关系^[21]；芦山地震与水库蓄水的关联度不大^[22]，而卡里巴水库的蓄水则导致了6.1级地震的发生^[23]。Albaric等^[24]研究南澳Paralana地区大规模水压致裂过程与地震释放的关系。Zang等^[25]研究欧洲地区地热工程对地震危险度的影响；Grogoli等^[26]和Kim等^[27]则指出，2017年11月韩国5.5级地震则受到地热系统2 a持续高压注水的影响，人类工业活动可能是造成该地震的重要因素；废水的注入和油气的开采也可能诱发地震的发生^[28]，美国Okalahoma地区油气开采注水可能造成中小地震数量增多^[29]。在这些问题中，数值模拟是解决问题的重要手段。
黏弹性有限元模型在当今地球动力学模拟中受到越来越多的重视，在各种地球动力学过程中得到了广泛应用。Freed和Lin^[30-32]利用黏弹性有限元模型研究地震同震和震后应力变化，认为1971年San FernandoM_W6.7地震可能触发了23 a后1994年Northridge M_W6.7地震；Wiseman等^[33]利用非均匀黏弹性模型研究2004年苏门答腊地震由于黏弹性松弛引起的震后变形。沈正康等^[34]采用黏弹性分层模型研究东昆仑活动断裂带大地震之间的黏弹性触发作用，认为中下地壳的黏性松弛使库仑应力变化随时间有增强的趋势。Hu等^[35]建立三维黏弹性有限元模型来研究1960年智利大地震的同震及震后变形，模拟结果较好地解释了大地震后35 a的GPS观测结果。Luo和Liu^[36]建立西藏东部及其邻区的三维黏弹塑性有限元模型，考虑重力、构造加载、下地壳和上地幔应力松弛、断层弱化等因素的影响，并模拟该区域历史大地震序列的时空演化，计算2008年汶川大地震的同震变形及同震库仑应力变化，探讨2008年汶川大地震的动力学成因。不同****利用黏弹性力学模型研究了冰川、地表水、地壳变形和地震之间存在复杂的相互作用，定量分析冰川载荷变化对地表变形和古今地震活动性变化的影响^[37-41]。
本文利用自主开发的Maxwell黏弹性有限元程序，研究1956—2106年间，在冰川不断消融引起的地表载荷和两侧地块的不断构造挤压作用下，祁连山地区地表变形、地应力的时空演化，定量评估研究区域主要断层带的地震危险性的时空变化。
1 地质背景与有限元模型为了研究由于冰川的消融及近南北向挤压对祁连山地区地表变形、地应力及主要断层带地震危险性变化的影响，我们建立了祁连山地区一条典型地质剖面的二维平面应变Maxwell黏弹性有限元模型^[42-43]。剖面左端点为(96°30′E, 37°N)，右端点为(98°30′E，39°30′N)，对应于图 1(a)中红色实线位置。
Fig. 1

	Download: JPG larger image
图 1 祁连山地质背景与剖面地形图[46] Fig. 1 Geological background and topographic section map of Qilian Mountains 图 1 祁连山地质背景与剖面地形图[46] Fig. 1 Geological background and topographic section map of Qilian Mountains -->

根据地震层析成像结果(CRUST1.0)建立有限元模型如图 2(a)所示，模型主要由3层介质组成，地表是根据实际地形插值绘制^[44]。莫霍面起伏参考地震层析成像的结果(CRUST1.0)，莫霍面以上2层为各向同性弹性介质，莫霍面以下为各向同性Maxwell黏弹性介质。为了突出断层带介质的差异，地质剖面上的5条断层看作宽度为300 m的断层带，断层的位置和几何产状来自前人的地质考察结果^[44]，视为横观各向同性介质^[45]。断层带介质比围岩要小一些，断层带介质的G₂为周围介质的一半。介质材料参数见表 1^[46]。本文的有限元模型采用四边形网格剖分，共有36 096个节点，35 700个单元，网格剖分见图 2(b)。
Fig. 2

	Download: JPG larger image
图 2 有限元模型与网格剖分图 Fig. 2 Finite element model and mesh model 图 2 有限元模型与网格剖分图 Fig. 2 Finite element model and mesh model -->

Table 1

表 1 黏弹性有限元模型的参数表[46]Table 1 Material parameters of the viscoelastic finite element model 材料 杨氏模量/Pa 泊松比ν 深度范围/km 黏度/ (Pa·s) 断层带 5.6e10 0.25 0~20 - 上地壳 5.6e10 0.25 0~20 - 下地壳 5.6e10 0.25 20~50 - 岩石圈地幔 7.1e10 0.25 50~100 1e20

表 1 黏弹性有限元模型的参数表[46]Table 1 Material parameters of the viscoelastic finite element model

边界条件：通过GPS插值得出剖面两侧的GPS速度值，将右侧边界设为水平方向位移分量为零，垂直方向剪应力为零；左侧边界水平速度为左右两侧GPS水平速度之差，垂直方向剪应力为零；下边界水平方向剪应力为零，垂直方向位移分量为零；上边界施加随时间变化的冰川融化卸载的载荷，相当于在上边界施加了向上的面力，其大小随冰川空间部位和融化时间而变化。值得注意的是，1956年祁连山地区的冰川范围比图 1(a)所示的范围要大很多。根据相关文献，本文假设冰川载荷的厚度和分布均随时间不断变化：冰川最低海拔随时间按10 m/a的速度上升，冰川厚度在4 500 m处的消融速率为0.32 m/a^[47]，从最高点到冰川最低海拔按线性分布，满足以下公式

$H = {H_0} \cdot \frac{{h - {h_{l\min }}}}{{{h_0} - {h_{l\min }}}}, $

(1)

其中：H为当前海拔冰川厚度，H₀为海拔4 500 m(h₀)处冰川厚度，h为当前海拔高度，h_lmin为冰川最低海拔。
计算时间阶段：1956—2106年，共150 a，计算的时间步长为1 a。由于Maxwell黏弹性模型是线性的，本文只考虑边界构造加载和地表冰川载荷变化的影响，不考虑初始应力场的影响，因此在模型中不考虑重力的因素。
2 模拟计算结果通过与1956年的初始位移场相减，得到1956—2106年间共150 a的位移应力变化量。图 3给出第50、第100和第150年位移的时空变化。1956—2106年间，随着时间推移，水平位移不断增加，主要受到模型两侧挤压变形的影响，冰川载荷的变化对其影响相对较小(图 3(a))；冰川载荷的影响主要表现在垂直位移上，在冰川载荷的逐渐消融和两侧挤压的共同作用下，地表呈现明显的向上抬升的现象，最大抬升量达0.55 m以上，主要位于海拔超过4 500 m的高海拔地区，海拔越高的地区，冰川载荷减小越多，抬升幅度越大，模型北段的抬升量不够明显(图 3(b))。
Fig. 3

	Download: JPG larger image
图 3 1956—2106年间位移变化等值线图(等值线单位：m) Fig. 3 Displacement contour map during 1956-2106 图 3 1956—2106年间位移变化等值线图(等值线单位：m) Fig. 3 Displacement contour map during 1956-2106 -->

1956—2106年间，水平应力、垂直应力和剪应力都发生了显著变化(图 4)。对于水平应力和垂直应力的变化，大于零代表拉应力增加，压应力减小，对应暖色区(红色)；小于零代表拉应力减小，压应力增加，对应冷色区(蓝色)，以下同。
Fig. 4

	Download: JPG larger image
图 4 1956—2106年间应力变化等值线图(等值线单位：Pa) Fig. 4 Horizontal stress contour map during 1956-2106 图 4 1956—2106年间应力变化等值线图(等值线单位：Pa) Fig. 4 Horizontal stress contour map during 1956-2106 -->

在构造挤压和冰川卸载的双重作用下，各主要断层带的正应力和剪应力呈现不同的演化趋势，最终会对地震危险性的改变带来不同影响，可以通过库仑应力变化进行展示。库仑应力变化的计算公式如下

$\Delta {\rm{CFS}} = \Delta + \mu \Delta {\sigma _n}, $

(2)

其中：Δτ为断层上剪应力变化量，以有利于地震发生方向为正；Δσ_n为断层正应力变化量，以拉应力为正，μ为摩擦系数^[48-49]。
图 5给出各主要断层带中间深度点处库仑应力变化ΔCFS在1956—2106年间随时间的变化曲线，可以清楚看出位于地质剖面不同部位的不同产状的各条断层，其库仑应力变化ΔCFS随时间的变化曲线差异性很大，也意味着不同断层带的地震危险性随时间的变化的差异性很大。祁连山北段断裂、哈拉湖南山断裂和哈拉湖盆地北缘断裂的库仑应力变化ΔCFS在这150 a期间都是大于零的，并且持续增加，地震危险性也持续增加，需要重点关注。右侧2条断层、特别是祁连山北断层ΔCFS基本变化不大，增长缓慢，非常接近零，可能反映仅仅在构造应力作用下，地震危险性增加的速率相对较小，而冰川的融化，特别在冰川融化量最大的最近50 a前后，对融化地区下方的断层地震危险性有较大影响。
Fig. 5

	Download: JPG larger image
图 5 各主要断层带中间深度点ΔCFS变化曲线 Fig. 5 Curves of ΔCFS with time at the center of major faults 图 5 各主要断层带中间深度点ΔCFS变化曲线 Fig. 5 Curves of ΔCFS with time at the center of major faults -->

3 讨论本文研究祁连山地区现代冰川融化对地表变形、地应力和断层地震危险性时空变化的影响。现代冰川的分布如图 1所示，文献表明其正在加速融化^[3]。根据计算结果，绘制模拟计算50、100、150 a和现今地表抬升速度图，如图 6(a)所示。图像表明，当地表存在冰川时，海拔越高，冰川融化速度越快，地表抬升速度越大。冰川融化速度越快的地区，地表抬升速度越大；当冰川融化后，地表抬升速度依赖于边界构造的挤压作用。2006年祁连山地区研究剖面的垂直抬升速率为1~5 mm/a，2009年抬升速率略有减小，这些模拟结果与前人观测得到祁连山地区垂直抬升速率^[50-51]较为一致。2056和2106年，冰川完全融化，抬升速率趋于稳定，最大值只有1 mm/a左右。
Fig. 6

	Download: JPG larger image
图 6 地表抬升速率与地震数据分析图 Fig. 6 Analysis of surface uplift rates and seismic data 图 6 地表抬升速率与地震数据分析图 Fig. 6 Analysis of surface uplift rates and seismic data -->

然后对研究区域的地震震级完备性进行验证。结果表明，图 1(a)所示的研究区域内1970—2016年间3.0~3.9级地震的数目为437，4.0~4.9级地震数目为254，5.0~5.9级地震数目为67，6.0~6.9级地震数目为13，7.0~7.9级地震数目为1，3级以上的地震基本满足古登堡和里克特地震震级频度G-R关系，3.5级以上的地震完备性较好，如图 6(b)所示。图 6(c)给出本文研究的祁连山剖面两侧30 km以内的1956—2016年的3级以上历史地震分布，这60 a的历史地震每隔10 a以不同颜色的圆圈表示。从空间上看，在剖面的南段，历史地震分布最多，主要位于剖面左一和左二的大柴旦—宗务隆山、哈拉湖南山2条断层之间，震源深度主要位于5~30 km。研究剖面的中段地震极少，北段有零星的地震分布。这种空间分布是地质构造条件决定的，不是本文讨论的范围。我们关心的是从时间来看，1996年以后的地震活动有显著增强，虽然不能肯定它们就是冰川融化造成的应力变化导致的，但这种地震活动性的时间变化与本文计算得到剖面左侧3条断层地震危险性增加的结果的确是较为一致的。
现代冰川的消融，导致本文研究剖面除左侧第1条断层之外的4条断层(从南到北依次为：HSF，哈拉湖南山断裂；HNF，哈拉湖盆地北缘断裂；YAF，阴凹槽断裂；QNF，祁连山北缘断裂)的断层面正应力增加(正应力以拉为正)，不仅有利于这4条断层逆断层地震的触发，而且有利于这4条断层发生左旋走滑断层地震的触发。祁连山地区的逆断层地震和左旋走滑地震都可能因为现代冰川融化而加速触发。
为了解现代冰川融化对祁连山地区地震危险性变化影响的程度，不妨比较一下历史大地震对祁连山地区主要断裂造成的库仑应力大小。2008年在青藏高原东缘和四川盆地交界的龙门山断层发生M_W7.9汶川大地震，汶川地震震中离祁连山地区的距离有800 km以上，造成祁连山地区主要断层的库仑应力变化接近于零，不易直接触发地震^[52]。张贝等^[53]利用全球分层介质研究2015年4月25日M_S8.1尼泊尔大地震的同震应力触发，对1 000 km以外的祁连山地区的应力变化影响很小。任天翔等^[54]利用三维含地形高程的横向不均匀性椭球型地球有限元模型研究2001年昆仑山口西M_S8.1大地震对周围断层的同震库仑应力变化的影响，结果表明祁连山地区的祁连山断裂的平均库仑应力变化为-0.001 kPa，宗务隆—南祁连冲断裂的平均库仑应力变化为-1.29 kPa，均不利于触发地震。通过这些比较，可以认识到本文计算的冰川融化与构造应力叠加对祁连山主要断层造成的库仑应力变化速率，最大可达+1 kPa/a以上，而且在150 a内不断累积，这是相当可观的效应，不容忽视。
4 结论本文利用自主开发有限元程序，构建祁连山地区的二维平面应变黏弹性有限元模型，几何模型考虑了实际地形和莫霍面起伏，模拟祁连山地区剖面上边界构造加载和冰川载荷不断消融对地表变形、地应力和断层地震危险性的时空变化，得出以下初步认识：
1) 在构造加载和冰川持续卸载的情况下，冰后回弹效应明显，最大垂直抬升可达0.55 m以上，海拔越高，位移回弹量也越大。随冰川载荷的减少，地表抬升幅度逐渐增大，海拔较高的中段抬升量最大。
2) 1956—2106年间，研究区域的水平应力、垂直应力及剪应力在冰川载荷卸载和边界构造加载的双重作用发生显著变化；各主要断层带的地震危险性呈现差异性的变化趋势，其中大柴旦—宗务隆山断层、哈拉湖南山和哈拉湖盆地北缘断层ΔCFS大于零，且随时间一直增大，地震危险性持续增加，应该引起足够重视；阴凹槽、祁连山北缘2条断层ΔCFS变化不大，未来发生地震的风险相对较小。
本文只考虑了二维平面应变黏弹性模型，未考虑祁连山地区现代冰川的三维空间展布，对冰川载荷的时空变化估计还比较粗略，下一步将开展祁连山地区三维黏弹性有限元模型的工作，以期进一步研究冰川消融对该地区地震危险性时空变化的影响。
感谢中国地震局地质研究所张竹琪博士在地震震源数据和绘图方面提供的帮助。
参考文献

[1]	Dobrowolski A B. Historia naturalna lodu[M]. Warszawa: Wyd kasa im Mianowskiego, 1923.
[2]	Cossart E, Mercier D, Decaulne A, et al. Impacts of post-glacial rebound on landslide spatial distribution at a regional scale in northern Iceland (Skagafj?reur)[J]. Earth Surface Processes and Landforms, 2014, 39(3): 336-350. Doi:10.1002/esp.3450
[3]	Sun M P, Liu S Y, Yao X J, et al. Glacier changes in the Qilian Mountains in the past half-century: based on the revised First and Second Chinese Glacier Inventory[J]. Journal of Geographical Sciences, 2019, 28(2): 206-220.
[4]	Yao T D, Thompson T, Yang W, et al. Different glacier status with atmospheric circulations in Tibetan Plateau and surroundings[J]. Nature Climate Change, 2012, 2: 663-667. Doi:10.1038/nclimate1580
[5]	Ding B H, Yang K, Yang W, et al. Development of a water and enthalpy budger-based glacier mass balance model (WEB-GM) and its preliminary validation[J]. Water Resources Research, 2017, 53: 1-33. Doi:10.1002/wrcr.22248
[6]	杨威, 周尚哲, Chiyuki N, 等. 祁连山中段新冰期以来的冰川序列及环境变化[J]. 兰州大学学报(自然科学版), 2006, 42(6): 12-15. Doi:10.3321/j.issn:0455-2059.2006.06.003
[7]	Zhu M H, Yao T D, Yang W, et al. Energy- and mass-balance comparison between Zhadang and Parlung No.4 glaciers on the Tibetan Plateau[J]. Journal of Glaciology, 2015, 61(227): 595-607. Doi:10.3189/2015JoG14J206
[8]	Yu W S, Yao T D, Kang S C, et al. Different region climate regimes and topography affect the changes in area and mass balance of glaciers on the north and south slopes of the same glacierized massif (the West Nyainqentanglha Range, Tibetan Plateau)[J]. Journal of Hydrology, 2013, 495: 64-73. Doi:10.1016/j.jhydrol.2013.04.034
[9]	葛伟鹏, 王敏, 沈正康, 等. 柴达木-祁连山地块内部震间上地壳块体运动特征与形变模式研究[J]. 地球物理学报, 2013, 56(9): 2994-3010.
[10]	刘杰. 基于GPS和水准数据的微粒群算法联合反演祁连山北缘断裂三维滑动速率[J]. 大地测量与地球动力学, 2019, 39(9): 906-909.
[11]	姜佳佳, 张辉. 祁连山地震带强震前ML4.0地震活动特征研究[J]. 数字地震观测与预报新方法及新技术, 2018, 2(8): 79-80.
[12]	姜佳佳, 张辉, 窦喜英, 等. 祁连山地震带及阿尔金地区MS ≥ 4地震活动增强研究[J]. 国际地震动态, 2019, 2(8): 72-73. Doi:10.3969/j.issn.0253-4975.2019.08.060
[13]	刘小凤, 梅秀苹, 冯建刚, 等. 青藏高原北部地区地震基本活动状态定量评价[J]. 西北地震学报, 2011, 33(2): 130-136.
[14]	刑武成, 李忠勤, 张慧, 等. 1959年来中国天山冰川资源时空变化[J]. 地理学报, 2017, 72(9): 1594-1605.
[15]	Liu C, Linde A T, Sacks I S. Slow earthquakes triggered by typhoons[J]. Nature, 2009, 459: 833-836. Doi:10.1038/nature08042
[16]	Bettinelli P, Avouac J P, Flouzat M, et al. Seasonal variations of seismicity and geodetic strain in the Himalaya induced by surface hydrology[J]. Earth and Planet Science Letters, 2008, 266: 332-344. Doi:10.1016/j.epsl.2007.11.021
[17]	Bollinger L, Perrier F, Avouac J P, et al. Seasonal modulation of seismicity in the Himalaya of Nepal[J]. Geophysical Research Letters, 2007, 34: L08304.
[18]	庞亚瑾, 张怀, 程惠红, 等. 华北地区地下水开采对地壳应力的影响[J]. 地球物理学报, 2016, 59(4): 1394-1402.
[19]	程惠红, 郑亮, 张贝, 等. 评"紫坪铺水库造成孔隙弹性耦合变化及其对2008年汶川地震触发作用"二维模拟的局限性[J]. 地球物理学报, 2015, 58(4): 1458-1461.
[20]	程惠红, 张怀, 石耀霖. 关于紫坪铺水库蓄水是否与汶川地震有关的影响因素综合分析[J]. 地球物理学报, 2015, 58(4): 1220-1235.
[21]	孙玉军, 张怀, 董树文, 等. 利用三维孔隙弹性模型探讨紫坪铺水库对汶川地震的影响[J]. 地球物理学报, 2012, 55(7): 2253-2361.
[22]	程惠红, 张怀, 石耀霖. MS7.0芦山地震的发生与水库有关吗?[J]. 科技导报, 2013, 31(31): 31-37. Doi:10.3981/j.issn.1000-7857.2013.31.004
[23]	程惠红, 张怀, 朱伯靖, 等. 卡里巴水库蓄水引起库区应力场变化影响分析[J]. 地震, 2013, 33(4): 32-42.
[24]	Albaric J, Oye V, Langet N, et al. Monitoring of induced seismicity during the first geothermal reservoir stimulation at Paralana. Australia[J]. Geothermics, 2014, 52: 120-131. Doi:10.1016/j.geothermics.2013.10.013
[25]	Zang A, Oye V, Jousset P, et al. Analysis of induced seismicity in geothermal reservoirs: an overview[J]. Geothermics, 2014, 52: 6-21. Doi:10.1016/j.geothermics.2014.06.005
[26]	Grigoli F, Cesca S, Rinaldi A P, et al. The novermber 2017 MW5.5 Pohang earthquake: a possible case of induced seismicity in South Korea[J]. Science, 2018, 360(6392): 1003-1006. Doi:10.1126/science.aat2010
[27]	Kim K H, Ree J H, Kim Y H, et al. Assessing whether the 2017 MW5.4 Pohang earthquake in South Korea was an induced event[J]. Science, 2018, 360(6392): 1007-1009. Doi:10.1126/science.aat6081
[28]	Keranen K M, Weingarten M, Abers G A, et al. Sharp increase in central Oklahoma seismicity since 2008 induced by massive wastewater injection[J]. Science, 2014, 345: 448-451. Doi:10.1126/science.1255802
[29]	Walsh F R, Zoback M D. Oklahoma's recent earthquakes and saltwater disposal[J]. Research article, 2015, e1500195: 1-9.
[30]	Freed A M, Lin J. Time-dependent changes in failure stress following thrust earthquakes[J]. Journal of Geophysical Research, 1998, 103: 24393-24409. Doi:10.1029/98JB01764
[31]	Freed A M, Lin J. Delayed triggering of the 1999 Hector Mine earthquake by viscoelastic stress transfer[J]. Nature, 2001, 411: 180-183. Doi:10.1038/35075548
[32]	Freed A M, Lin J. Accelerated stress buildup on the southern San Andreas fault and surrounding regions caused by Mojave Desert earthquakes[J]. Geology, 2002, 30: 571-574. Doi:10.1130/0091-7613(2002)030<0571:ASBOTS>2.0.CO;2
[33]	Wiseman K, Burgmann R, Freed A M. Viscoelastic relaxation in a heterogeneous Earth following the 2004 Sumatra-Andaman earthquake[J]. Earth and Planetary Science Letters, 2015, 431: 308-317. Doi:10.1016/j.epsl.2015.09.024
[34]	沈正康, 万永革, 甘卫军, 等. 华北强震序列中积累应力变化场演化及地震发生研究[C]//中国地球物理学会第十九届年会论文集. 南京, 2003: 276.
[35]	Hu Y, Wang K, He J, Klotz J, et al. Three-dimensional viscoelastic finite element model for postseismic deformation of the great 1960 Chile earthquake[J]. Journal of Geophysical Research, 2004, 109: B12403. Doi:10.1029/2004JB003163
[36]	Luo G, Liu M. Stress evolution and fault interactions before and after the 2008 Great Wenchuan earthquake[J]. Tectonophysics, 2010, 491: 127-140. Doi:10.1016/j.tecto.2009.12.019
[37]	Peltier W R. Ice age paleotopography[J]. Science, 1994, 265: 195-201. Doi:10.1126/science.265.5169.195
[38]	Stewart I S, Sauber J, Rose J. Glacio-seismotectonics: ice sheets, crustal deformation and seismicity[J]. Quaternary Science Reviews, 2000, 19: 1367-1389. Doi:10.1016/S0277-3791(00)00094-9
[39]	Wahr J, Velicogna I. What might GRACE contribute to studies of post glacial rebound[J]. Space Science Reviews, 2003, 108(1): 319-330. Doi:10.1023/A:1026183526762
[40]	Wang H S, Jia L L, Steffen H, et al. Increased water storaage in North America and Scandinavia from GRACE gravity data[J]. Nature Geoscience, 2013, 6(1): 38-42. Doi:10.1038/ngeo1652
[41]	Dietrich R, Ivins E R, Casassa G, et al. Rapid crustal uplift in Patagonia due to enhanced ice loss[J]. Earth and Planetary Science Letters, 2010, 289(1/2): 22-19.
[42]	Zheng G, Wang H, Wright T J, et al. Crustal Deformation in the India-Eurasia Collision Zone From 25 Years of GPS Measurements[J]. Journal of Geophysical Research: Solid Earth, 2017, 122(11): 9290-9312. Doi:10.1002/2017JB014465
[43]	李虹蓉, 杨太保, 田洪阵, 等. 1987-2015年祁连山冰川退缩及未来变化趋势研究[J]. 安徽农业科学, 2017, 45(30): 59-63. Doi:10.3969/j.issn.0517-6611.2017.30.019
[44]	Zhang H P, Zhang P Z, Prush V, et al. Tectonic geomorphology of the Qilian Shan in the northeastern Tibetan Plateau: insights into the plateau formation processes[J]. Tectonophysics, 2017, 706-707: 103-115. Doi:10.1016/j.tecto.2017.04.016
[45]	Hu C B, Cai Y E, Wang Z M. Effects of large historical earthquakes, viscous relaxation, and tectonic loading on the 2008 Wenchuan earthquake[J]. Journal of Geophysical Research, 2012, 117: B06410.
[46]	Grollimund B, Zoback M D. Post glacial lithospheric flexure and induced stresses and pore pressure changes in the northern North Sea[J]. Tectonophysics, 2000, 327: 61-81. Doi:10.1016/S0040-1951(00)00162-1
[47]	Gardner A S, Moholdt G, Cogley J G, et al. A reconciled estimate of glacier contributions to sea level rise: 2003 to 2009[J]. Science, 2012, 340: 852-857.
[48]	Oppenheimer D H, Reasenberg P A, Simpson R W. Fault plane solution for the 1984 Morgan Hill, California, earthquake sequence: Evidence for the state of stress on the Calaveras Fault[J]. Journal of geophysical research, 1988, 93(B8): 9007-9026. Doi:10.1029/JB093iB08p09007
[49]	Harris R A. Introduction to special section: Stress triggers, stress shadows, and Implications for seismic hazard[J]. Journal of Geophysical Research, 1998, 103(B10): 25, 347-24, 358.
[50]	刘峡, 黄立人, 杨国华, 等. 北祁连-河西走廊地区垂直形变与构造应力场关系的初步研究: 三维有限元拟合结果分析[J]. 地震地质, 2003, 25(2): 307-316. Doi:10.3969/j.issn.0253-4967.2003.02.016
[51]	梁诗明. 基于GPS观测的青藏高原现今三维地壳运动研究[D]. 北京: 中国地震局地质研究所, 2014.
[52]	Luo G, Liu M. Stressing rates and seismicity on the major faults in Eastern Tibetan Plateau[J]. Journal of Geophysical Research: Solid Earth, 2018, 123(12): 10968-10986. Doi:10.1029/2018JB015532
[53]	张贝, 程惠红, 石耀霖. 2015年4月25日尼泊尔MS8.1大地震的同震效应[J]. 地球物理学报, 2015, 58(5): 1794-1803.
[54]	任天翔, 程惠红, 张贝, 等. 2001年昆仑山口西MS8.1地震对周围断层的应力影响数值分析[J]. 地球物理学报, 2018, 61(12): 4938-4850.