冰冻圈是指地球上存在的所有冰体所构成的空间，包括包含冰的云层、地下冰、积雪、水冰、冰川、冰盖等瞬时、短期、季节和永久存在的冰体，最早在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)中红色实线位置。

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

边界条件：通过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))。

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

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

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

(2)

其中：Δτ为断层上剪应力变化量，以有利于地震发生方向为正；Δσ_n为断层正应力变化量，以拉应力为正，μ为摩擦系数^[48-49]。

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

3 讨论

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

然后对研究区域的地震震级完备性进行验证。结果表明，图 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变化不大，未来发生地震的风险相对较小。

本文只考虑了二维平面应变黏弹性模型，未考虑祁连山地区现代冰川的三维空间展布，对冰川载荷的时空变化估计还比较粗略，下一步将开展祁连山地区三维黏弹性有限元模型的工作，以期进一步研究冰川消融对该地区地震危险性时空变化的影响。

感谢中国地震局地质研究所张竹琪博士在地震震源数据和绘图方面提供的帮助。