0 引言

冰冻圈是指地球上所有自然存在的冰体所形成的圈层，包括瞬时、短期、季节和永久存在的冰体(谢自楚和刘潮海，2010).冰川和冰盖作为冰冻圈的主体，占陆地总面积的11%，在调节地球气候方面，起到了至关重要的作用.冰后回弹和冰川均衡调整对于地球各圈层产生重要影响，可能引起地形地貌的巨大变化，甚至引起地球扁率的变化(Kaufmann and Lambeck, 2000；Dickey et al., 2003；Bourgois et al., 2016；de Quay et al., 2019；孙付平，1995；闫昊明等，2004；谢自楚和刘潮海，2010)，可能引发大地震(Spada et al., 1992；Arvidsson，1996；Wu et al., 1999；Bungum et al., 2010；Fjeldskaar et al., 2000)，也可能引发滑坡等地质灾害(Cossart et al., 2014).冰川载荷的短期和长期地质效应正在受到越来越多的关注.

随着全球气温的不断升高，地球上的冰川正在加速消融.从整个地质历史来看，冰川增厚和消融是交替出现的，冰川载荷随地质时间不断改变.由于冰川载荷的时空变化，固体地球内部发生冰川均衡调整现象，并引起地表的冰后回弹.冰川均衡调整(Glacial Isostatic Adjustment，GIA)主要是由地球的冰川载荷和海水负荷的变化引起，表现为地幔物质的流动和地壳的变形(汪汉胜等，2009).冰后回弹(Postglacial Rebound)主要是指冰川消融、冰川载荷减小甚至消失所引起的固体地球的即时弹性变形和冰后黏弹性变形，最终造成地表的随时间变化的不均匀抬升或沉降现象(崔立鲁和朱贵发，2015).

地球介质的流变性模型对研究冰川消融引起的冰后回弹影响极大.空间大地测量技术的最近观测发现实际的回弹位移比理论模型预估值要大得多(朱心慧和孙付平，2005).卫星技术的应用可以发现冰后回弹的更多细节，发现不仅存在地表垂直速度的变化，也存在地表水平速度的变化(James and Lambert, 1993；James et al., 2000；Wahr and Velicogna, 2003).由于地球介质复杂的特殊的流变性，冰后回弹不仅会造成垂向应力和位移的变化，也会造成水平应力和位移的变化(孙付平等，1997；于男等，2014).因此，需要发展更符合实际的地球介质模型，考虑地球介质流变结构的不均匀性.

不同学者对地球典型的大型冰盖或冰川消融引起的冰后回弹做了大量的研究工作.如南极冰盖季节性降雨的质量平衡、冰盖消融引起的冰后回弹、海平面上升等问题都做了系统的研究(翟宁等，2009；崔立鲁和朱贵发，2015；史红岭等，2015；Shepherd, et al., 2012)，这些结果充分利用了卫星高程测量、干涉测量和重力测量的数据以及地球介质的流变性模型.格陵兰岛和南极冰盖边缘的冰川正在加速融化，它们对全球海平面上升的贡献正在越来越大.人们越来越希望能够获取冰川融化的真实数据(Pritchard et al., 2009).阿拉斯加东南部区域现今的地壳变形和地表垂直位移变化可能由约200年前的小冰期之后的冰川消退所控制和影响(刘静等，2018；Larsen et al., 2005).近几百年以来，全球气候变暖引起的冰川消退现象也在影响着南美洲西部的区域隆升和地幔对流(刘静等，2018).

Haskell早在1935年就提出了一个关于冰后回弹的解析解：一个黏性圆柱体在施加的垂向载荷移除后表面随时间逐渐抬升的位移解析解，可以用于解释固体地球在冰川消融后地表逐渐抬升的现象(Haskell，1935).Turcotte和Schubert用Matlab程序进行积分运算，绘制了Haskell(1935)的解析解曲线，对冰后回弹进行了研究(Turcotte and Schubert, 2002).Peltier通过对应原理，推导出了Maxwell黏弹性体的地球在地表块状载荷作用下地表抬升或沉降的解析解(Peltier，1974).这些解析解较为简单，是认识冰后回弹的基本动力学特征的研究基础.

目前，数值模拟方法已经成为研究冰后回弹和冰川均衡调整的重要手段.冰后回弹的时间尺度既不属于完全弹性尺度，又不同于流体的长时间尺度，因此，对于冰后回弹问题来说，黏弹性模型是十分合适的(肖强和许厚泽，1990).不同学者利用黏弹性力学模型研究了冰川、地表水、地壳变形和地震之间存在复杂的相互作用(Peltier，1994；Stewart et al., 2000；Wahr and Velicogna, 2003；Wang et al., 2012；Dietrich et al., 2010).张怀等利用自主开发的一套三维并行冰盖演化模型模拟了格陵兰岛冰盖的演化过程，考虑了复杂地形起伏和不均匀冰盖厚度、全球气候变化对冰盖演化的影响(Zhang et al., 2011)；冷伟等利用三维非线性冰川演化有限元模型模拟冰川应力场和速度场的演变(Leng et al., 2012).

北欧是研究冰后回弹和冰川均衡调整的重要区域.不同研究者积累了北欧地区第四纪冰后期以来的冰盖厚度和范围的记录，气温变化、平均海平面变化、古地震时空分布、现今地震活动性分布、现今地震震源机制解、断层分布，还积累了大量的地应力测量的大量资料，为深入研究冰后回弹和冰川均衡调整的动力学过程，尤其是地球介质的流变性特征、地应力时空演化与古地震、现今地震的关系提供了较为充分的资料.

Patton等(2017)建立了末次大冰期前的欧洲冰盖的三维热力耦的斯托克斯流体模型，包括初、边值条件，空间网格分辨率为10 km，研究了Celtic、Fennoscandian、Barents Sea三个冰盖自37~19 ka以来冰盖扩展加厚的过程，以及19~9 ka以来的冰盖退冰过程.Wu等(2005)利用18层黏弹性有限元模型研究了冰后回弹效应，结果表明，陆地垂直速率、海平面垂直速率、水平速度、大地水准高度等的观测结果有助于探测地幔的横向不均匀性.

Wolf(1987)在首先假定软流圈黏度为1×10²¹ Pa·s的情况下，利用半无限空间三层不可压缩Maxwell黏弹性半解析解和冰后回弹观测资料，推断出芬诺斯坎迪亚地盾(Fennoscandia)岩石圈厚度的上限为80 km，岩石圈黏度为1.2×10¹⁹ Pa·s.Johnston等(1999)利用5层不可压缩黏弹性模型研究了水平应力变化与冰川载荷的关系.Klemann和Wolf (1998)研究了冰川载荷分布、岩石圈厚度、软流圈黏度等因素对芬诺斯坎迪亚地盾冰后回弹应力变化的影响，结果表明岩石圈的黏弹性松弛的影响非常大.Wu等(1999)利用黏弹性模型研究了芬诺斯坎迪亚地盾冰后回弹引起的应力变化与研究区域断层不稳定性之间的关系，认为冰后回弹是在冰后期观测到的大的逆断层的可能成因.这些研究者大部分使用的是解析解或半解析解进行研究.

关于北欧的末次冰期之后退冰过程，有一个现象引人关注：那就是虽然现今地震活动性较弱，震级小，频度低，但古地震活动性强，震级大，频度高(Arvidsson and Kulhánek, 1994；Fjeldskaar et al., 2000；Mörner, 2004, 2005).相伴的现象就是现今冰后回弹垂直速率仅为几个mm·a^-1，但历史上地表垂直速率最大时可达10~50 cm·a^-1(Ekman and Mäkinen, 1996; Mörner, 2004；Olesen et al., 2013).末次冰期以来北欧地区的冰盖厚度加厚或减薄交替出现(Stroeven et al., 2016).虽然有研究者试图对地震活动性变化作定性的解释(Arvidsson and Kulhanek, 1994)，但定量化的解释，特别是因为冰盖反复加载、卸载的复杂时空演化历史下的定量化解释，迄今尚未见到.本文利用自主开发的研究冰后回弹和冰川均衡调整的黏弹性有限元程序，在考虑重力和构造加载、地球介质弹性的纵向和横向不均匀性以及黏性的分层性条件下，研究了末次冰期以来冰盖载荷变化引起的斯堪的纳维亚地壳变形及应力状态.对挪威北海北部1.1 Ma以来的冰盖载荷变化、特别是两万多年前冰盖消退引起的地表冰后回弹的应力变化进行了数值模拟，模拟结果与现今实际观测有很好的一致性，从而可以进一步定量化解释古今地震活动性差异的转变机制.

1 挪威北海现今地震活动性和地表变形特征

挪威北海北部现今的地震活动性与GPS变形如图 1所示.图中地震选自FENCAT地震目录的1980-2011年的地震，可以看出在挪威北海北部的现今地震主要集中于两个区域，一个位于海岸线以西的海洋一侧，一个位于海岸线以东的大陆一侧.由震源机制解所指示的水平P轴(压缩)在挪威北海北部海岸线以西是NWW-SEE为主，在大陆内部以NEE-SWW，接近EW水平方向为主，逆断层地震较为发育，大陆一侧的浅部地震还存在较为明显的水平T轴(拉张)，有少量的正断层地震(Arvidsson and Kulhánek, 1994; Fjeldskaar et al., 2000；Keiding et al., 2015).现今GPS速度资料表明，斯堪的纳维亚半岛的地表水平速度呈现明显的放射状分布，即由半岛中心向外部运动.现今GPS垂直速度分布还表明，现今地表垂直速度也是以半岛中心最大，可达8 mm·a^-1，垂直速率向四周逐渐减小，在海岸线附近只有1~2 mm·a^-1(Fjeldskaar et al., 2000；Grollimund and Zoback, 2000；Keiding et al., 2015).利用世界应力地图的最新资料发现挪威北海北部的最大主压应力的分布在方位角80~100°之间变化，接近东西向分布.

芬诺斯坎迪亚地盾的1.3万以来的古地震分布非常活跃，有35个震级大于6级的古地震分布，其中有8级以上的大地震(Mörner, 2004)，而现今地震的最大震级一般不超过4.2~4.3级，还发生过5级左右的历史地震，如1497年(4.8)、1759年(5.3)，1904年(5.4)，现今地震的活动性、震级、频度都不及古地震(Mörner, 2004).

末次冰期以来冰盖载荷变化引起的斯堪的纳维亚半岛自末次冰期以来经历复杂的冰盖加载、卸载历史(Mangerud et al., 2016; Stroeven et al., 2016).本文选取的挪威北海北部的一条近东西向的二维剖面，约从北纬61.5°，西经0.2°到北纬60.2°东经10.2°(Grollimund and Zoback, 2000)，冰盖厚度自1.1 Ma以来随时间的变化见图 2.现今海岸线约在图 2中的水平距离280 km处.自1.1 Ma以来，挪威北海北部冰盖厚度随时间和空间有很大的变化，1.1 Ma~11万年以前，冰盖基本只存在于大陆一侧，冰盖最大厚度约为1200 m(图 2a)；11~2万年以前，冰盖厚度迅速增加，尤其是海洋一侧，冰盖厚度超过1500 m，而大陆一侧420 km附近的冰盖厚度有所降低，从1200 m降到800 m左右；自2万年来以来，海洋一侧的冰盖迅速消失，大陆一侧的冰盖厚度也逐渐减小，从2万年前的1200 m减小到1万年前的800 m，本文还考虑了最终冰盖完全消失的情形(图 2f).

2 挪威北海北部冰后回弹引起的地壳变形和应力变化的计算

北欧地区自2万年以来，冰盖载荷的时空变化较大，是研究冰后回弹和冰川均衡调整的重要区域.我们利用自主开发的Maxwell黏弹性程序研究挪威北海北部自1.1 Ma以来，尤其是2万年以来冰盖载荷演化对冰后回弹和冰川均衡调整的影响(Grollimund and Zoback, 2000).模型描述如下：模型选取的是挪威北海北部的一条近东西向的二维剖面，约从北纬61.5°，西经0.2°到北纬60.2°东经10.2°(Grollimund and Zoback, 2000)，水平长度600 km，深度120 km，由5层介质构成，其中上地壳两层为弹性结构，下地壳及上地幔之上部分为两层黏弹性介质，上地幔之下半部分为不可压缩的黏弹性介质.模型左右两侧和底边均为法向方向位移固定(垂向位移为零)，切向方向自由滑动(剪应力为零)，地表垂向受冰盖载荷正压力作用(见图 2)，切向自由(剪应力为零).介质材料参数见表 1.

本文采用二维平面应变Maxwell黏弹性模型，本构方程如下(许鹤华等，2002；胡才博，2009)：

(1)

其中{σ}为应力，为{σ}对时间的导数，为{ε}应变对时间的导数，K、G和η分别为体积模量、剪切模量和黏度.

通过比较数值解和黏弹性解析解的位移和应力结果，验证自主开发的二维黏弹性有限元程序的正确性.黏弹性问题解析解见付真(2008)；胡才博(2009)，位移和应力解析解如下：

(2)

(3)

其中ν为泊松比，

利用有限元程序进行计算，模拟计算结果与解析解理论公式结果进行对比，如图 4所示.

可以看出，无论是位移还是应力的数值模拟结果，自主开发的Maxwell黏弹性程序都具有不错的精度，和解析解非常吻合，为研究冰后回弹和冰川均衡调整奠定了坚实基础.

针对挪威北海北部建立的二维平面应变有限元模型(图 3)，我们先给模型施加重力让其达到重力均衡状态(图 5).这个初始场是从1.1 Ma之前计算到11万年之前，对应于图 2a的末尾阶段.为了描述的方便，本文都以离模型最左端(最西端)的水平距离和深度来指示位置.本文中，20 km深度以上是完全弹性介质，20 km以下是黏弹性介质.可以看出，重力均衡之后的水平正应力和垂直正应力都是成层分布的，并且两者在20 km以下趋于一致，而剪应力几乎为零，表现为静岩压力状态.我们在这个初始应力状态下，按照图 2所示的冰盖厚度随时间的变化，逐渐改变地表冰盖载荷，以后的每一步结果都可以通过与重力均衡的结果相减得到冰盖加载和卸载对地球介质产生的影响.

为了探讨挪威北海北部地震活动性变化的原因，我们绘制了11万年以来不同时间节点的水平应力、垂直应力和剪应力对比图(图 6、图 7和图 8).本文中，正应力以拉为正，正应力变化大于零，表示拉应力增大，压应力减小，对应于图 6和图 7中黄橙红等暖色区域；正应力变化小于零，表示拉应力减小，压应力增大，对应于图 6和图 7中青蓝等冷色区域.可以明显地看出，在深度35 km和50 km处，三组应力分量都存在较为明显的突变，反映了界面处材料突变和有限元网格突变的影响，其中水平应力在50 km深度处的突变最为显著.

挪威北海北部自11万年以来由于冰盖载荷的变化，水平应力随时间和空间的变化较为复杂，11~2万年之前，整条剖面除了水平距离420 km处的一段，主要对应于冰盖的加载阶段，相比于初始状态，20 km深度以上的弹性层，分为四个区域，其中水平压应力增加区位于水平距离220 km处，水平压应力减小区位于水平距离420 km处，且主要在10 km以浅的区域(暖色区域).剖面左端(西段)水平压应力减小，剖面右端(东段)水平压应力增大(图 6a，冷色区域).20~35 km深度的黏弹性层水平应力变化不大，而35~50 km深度的黏弹性层积累了较大的水平应力变化，幅值可达40 MPa，原因主要是因为后者的杨氏模量、密度和黏度比前者要大很多(见表 1).

11~2万年之间，冰盖载荷没有发生变化，但由于20 km之下的下地壳和上地幔的黏弹性松弛，到2万年之前，20~35 km的黏弹性层水平应力由于较大的变化，从而引起20 km的弹性层底部水平应力也有了很大的变化，其中水平距离220 km附近的一段水平压应力减小，420 km附近的一段水平压应力增加.10 km以浅，水平压应力增加区位于220 km处，水平压应力减小区位于420 km处，且主要在10 km以浅的区域，与之前的分布一致，但幅值有所减小(图 6b).

1.5万年之前，由于冰盖载荷持续的卸载，尤其是250 km以西的海洋一侧，冰盖完全消失，其余部分，除了420 km附近载荷增加，冰盖载荷大部分都是减小的，1.5万年之前的剖面20 km以浅的弹性层，水平应力的分布相较于11~2万年之前发生很大的变化，220 km的部分(海洋段)水平压应力减小，而420 km附近的部分(大陆段)水平压应力增大(图 6c).相比图 6a和b，35 km以下的黏弹性层水平应力得到释放，逐渐转移到20~35 km的黏弹性层，最终导致20 km以浅的弹性层的水平应力的分布发生极性变化(图 6c).

1.1万年之前，海洋一侧的载荷完全消失(280 km以西的部分)，大陆一侧(280 km以东)冰盖载荷继续卸载，水平压应力减小区比1.5万年之前继续扩大，覆盖380 km以西的大部分区域，400 km以东的大陆区域水平压应力增加(图 6d).1万年之前，因为时间相隔较小，冰盖载荷的卸载也不多，所以结果与1.1万年之前的比较相似，只是幅值有所减少(图 6e).当冰盖完全消融时，海洋一侧(280 km以西)，20 km以浅的弹性层水平压应力大部分都是减小的，但220 km附近的一段，3 km以浅水平压应力增加的；大陆一侧，380 km以东，大部分区域水平压应力仍然是增加的，只是增加的幅度比1万年之前要小一些，水平距离280~400 km在5 km以浅的部分区域，水平压应力减小，在浅部可能会发育小规模的正断层(图 6f)，这个结果与现今地应力测量显示的深部压缩为主、浅部存在拉张、断层分布、现今地震震源机制解的结果一致(Fjeldskaar et al., 2000；Arvidsson and Kulhánek, 1994; Arvidsson, 1996; Grollimund and Zoback, 2000; Olesen et al., 2013).

自2万年以来，随着冰盖载荷的持续卸载，水平距离400~500 km的弹性层底部水平压应力开始增加(图 6b)，并且这个水平应力增加的区域，随着时间演化，水平应力增加的区域逐渐扩大，从弹性层底部逐渐发展到地表，1.1万年前达到最大区域(图 6d)，然后水平应力增加的区域开始减小(图 6e)，但直到冰盖完全消融，依然存在水平应力增加的区域(图 6f)，这有利于大陆一侧的逆断层的发育和活化，这个模拟结果与末次冰后期古地震活动性较强、现今地震活动性较弱是一致的(Arvidsson and Kulhánek, 1994；Mörner, 2004, 2005; Olesen et al., 2013).

图 7给出了11万年以来垂直应力的变化.11~2万年，冰盖载荷绝大部分都是增加的，只有局部(420 km附近)是卸载的，这一期间冰盖载荷增加最多的部分在海洋部分(280 km以西)，因此垂直压应力在海洋部分是增加的(280 km以西)，增加的幅值可达10 MPa.420 km附近由于冰盖载荷减小，这一部分垂直压应力是减小的(图 7a，7b).自2万年以来，海洋部分冰盖载荷基本消失，陆地部分的冰盖载荷持续减小直至完全消融，垂直压应力在海洋部分基本不变，在大陆部分持续降低(图 7c, d, e, f)，大陆一侧垂直压应力的持续减小更有利于逆断层的发育和活化，这个结果也与现今地应力测量、断层分布、现今地震震源机制解的结果一致(Fjeldskaar et al., 2000；Arvidsson and Kulhánek, 1994; Arvidsson, 1996; Grollimund and Zoback, 2000; Olesen et al., 2013).

我们绘制了冰盖不同加载、卸载阶段剪应力τ_xy随时间、空间的演化图(图 8).可以看出，与水平正应力、垂直正应力变化(图 6，7)相比，因冰盖加载卸载引起的剪应力变化幅度较小，在5 MPa以内变化，说明主应力基本还是水平和垂直方向.当11万年前，冰盖加载时(图 2b初始)，剪应力是非均匀分布的，海岸线附近的大陆一侧，剪应力变化接近5 MPa，对原有的应力状态有加大改变(图 8a中的红色区域)，远离海岸线的海洋一侧和大陆一侧还存在两个方向相反的剪应力变化较大的区域(图 8a中的蓝色区域)，主要由地表冰盖载荷的非均匀分布而引起.随着下地壳和上地幔的黏弹性松弛，到2万年前(图 2c的末尾)，剪应力变化幅值有了很大的减小(图 8b).2万年以来，冰盖载荷基本以卸载为主，剪应力变化较大的区域主要集中在海岸线附近区域(图 8c-e)，其中1.5万年以前的海岸线附近的剪应力变化的方向发生了变化，由红色区域变成了蓝色区域(图 8c)，1.1万年以前和1万年以前，该区域的剪应力变化幅值逐渐增大，接近5 MPa(图 8d，e)，当冰盖完全消失时，该区域的剪应力变化幅值迅速减小至零左右.2万年以来的冰盖卸载所对应的海岸线附近的剪应力幅值先增大后减小，1.5万年以来的剪应力变化增加的阶段与古地震活动性增强、频度高、震级大相对应，而当冰盖完全消融时，剪应力变化基本降至零，与现今地震活动性较弱、频度低、震级小相对应(Arvidsson and Kulhánek, 1994；Mörner, 2004, 2005; Olesen et al., 2013).

我们的数值模拟得到的现今地表垂直变化速率见图 9a.1.5万年以来的地表垂直变化速率是在不同的冰盖载荷阶段是不同的，存在明显的差异性垂直位移变化，其中大陆一侧地表垂直变化速率在1~4 mm·a^-1，明显地高于海洋一侧(模型左侧，280 km以西)，大陆一侧的现今地表垂直变化速率与观测数据比较吻合(图 9a中的黑色曲线)(Grollimund and Zoback, 2000; Fjeldskaar et al., 2000).我们还给出了1.5万年以来该剖面水平变化速率随时间的变化图 9b)，水平变化速率较小，绝对值在2 mm·a^-1之下，其中1万年前的水平变化速率最大(图 9b中的蓝色曲线)，反映了当时地表冰川载荷剧烈变化的影响，现今该剖面水平变化速率最小，低于1 mm·a^-1(图 9b中的黑色曲线)，与该剖面现今冰川完全消失相对应，同时反映出当前该地区的构造区域作用不强.

根据图 6和图 7，20 km深度以上的弹性层自2万年以来受冰盖载荷变化和黏弹性层松弛的影响，形成了两个应力变化较为剧烈的区域，分别位于水平距离220 km和420 km两个部分，分别对应海洋段和大陆段.我们取深度10 km，水平距离分别为220 km和420 km两点来具体分析水平应力和垂直应力的变化(表 2).这两个点分别与挪威北海北部的两条地震带的位置对应.表 2中的负值表示压应力增加，正值表示压应力减小.对应海洋段的水平距离220 km处的点，11~2万年之前，该点之上的冰盖是急剧加载的，因此这段时间该点水平压应力是增加的，增加的幅度为11.03~3.38 MPa，垂直压应力也是增加的，增加的幅度为17.31~17.27 MPa.1.5万年之前，该点水平压应力减小11.06 MPa，1.1万年9.02 MPa，1.0万年减小10.66 MPa，至冰盖完全消融时减小8.74 MPa，而垂直压应力自1.5万年以来变化不大，变化幅度约0.1 MPa.

对应大陆段的420 km，水平压应力从11万年前减小8.17 MPa，2万年前增加2.00 MPa，1.5万年增加13.74 MPa，1.1万年增加12.76 MPa，1.0万年增加12.84 MPa，至冰盖完全消融时增加9.21 MPa，该点水平压应力存在一个明显的先增加后减小的过程，这与古地震活动性较强(1.5~1万年之前)，而现今地震活动性较弱是一致的(Arvidsson and Kulhánek, 1994；Mörner, 2004, 2005).该点的垂直压应力基本都是减小的，对应于大陆段的冰盖卸载，变化幅度在2.44~4.94MPa，至冰盖完全消融时，该点垂直压应力将减小至12.37 MPa，垂直压应力的减小有利于该处逆断层的生成或活化，这个结果也与现今地应力测量、断层分布、现今地震震源机制解的结果一致(Fjeldskaar et al., 2000；Arvidsson and Kulhánek, 1994; Arvidsson, 1996; Grollimund and Zoback, 2000).

相比正应力变化，剪应力变化较小，最大值为水平距离为420km处的大陆段10 km深度处，在1.1万年以前的剪应力变化幅度为2.8 MPa，之后减小到1 MPa以下(表 2).以上结果表明，自2万年以来，由于冰盖载荷的时空变化，这两个点的应力状态也跟随发生变化，导致古地震和现今地震的地震活动性存在显著差别.

3 讨论和结论

本文利用黏弹性模型研究了研究挪威北海北部末次冰期以来的冰后回弹和应力演化过程，主要考虑了两方面因素的影响.

一是末次冰期11万年以来复杂的冰盖加载、卸载的过程，其中海岸线以西的海洋部分在2万年以来基本上冰盖完全卸载，而海岸线以东的大陆部分，大部分区域都是卸载的，但局部区域存在加卸载反复变化的情形，增加了冰后回弹和应力演化过程的复杂性.前人利用解析解或半解析解研究冰盖加卸载对冰后回弹的影响，冰盖加卸载过程较为简单，只能给出一些定性的分析(Haskell，1935；Wolf，1987；Johnston et al., 1999；Wu et al., 1999).冰盖加卸载过程的复杂性，很难用解析解或半解析解来进行研究，这也是本文采用有限元数值模拟的原因.

二是我们使用的Maxwell黏弹性模型，可以考虑在复杂冰盖载荷变化过程中应力在不同的弹性层和黏弹性层之间的转移，这种应力转移体现了弹性、黏弹性介质的特性，黏弹性层的应力总体是松弛的，应力逐渐转移到上覆弹性层的底部并逐渐向浅部发展，造成弹性层深部和浅部应力状态的差异性分布.我们的模拟结果显示海岸线以东的大陆一侧自2万年以来，弹性层深部水平压应力增大，有利于逆断层的发育或活化，与前人的结果一致(Wu et al., 1999)，同时弹性层的浅部水平压应力还可能减小，有利于浅表正断层的发育或活化，这与现今地表抬升引起局部正断层的分布是一致的(Fjeldskaar et al., 2000；Grollimund and Zoback, 2000).如果全部采用弹性层来研究复杂冰盖载荷变化引起的冰后回弹和应力演化，则只能即时反映地表冰盖载荷的变化，而不能反映下地壳和上地幔的黏弹性松弛效应(Grollimund and Zoback, 2000).

在考虑了复杂的冰盖加卸载过程以及黏弹性流变结构的基础上，我们进行了有限元数值模拟，得到以下初步的结论.

(1) 自2万年以来的冰盖卸载，引起挪威北部北海海岸线附近的剪应力变化先增加至1万年以前的5 MPa的最大值，又迅速减小至现今冰盖完全消失时的零值附近，大陆一侧的水平压应力在1.1~1万年以前增大到最大值，垂直压应力持续降低，可以用来解释古地震和现今地震活动性的显著差异.

(2) 模拟得到海岸线附近的水平应力和垂直应力的两个应力集中区，与挪威北海北部海岸线附近的海洋和大陆的两条地震带是对应的.模拟得到的地表抬升速率在不同的冰盖载荷阶段有不同的分布，现今地表抬升速率与观测结果一致，现今抬升向上最快的地方，有利于浅表次要的正断层活动，对应于模拟得到的上地壳浅表的水平应力拉张区；

(3) 我们自主开发一套Maxwell黏弹性有限元程序，并通过了解析解的验证，可以很好地模拟不同冰盖载荷情况下冰后回弹和冰川均衡调整的科学问题.目前，我们的模型还较为简单，没有考虑冰盖消融所引起的海平面上升的影响，没有考虑更复杂地球流变模型(如双黏度Burgers模型)的影响，没有考虑三维冰盖消融的影响，今后将在这些方面做更深入细致的工作.

致谢  本研究得到北京大学蔡永恩教授、中国科学院大学皇甫鹏鹏博士很多的意见和建议，两位匿名评审人对初稿和二审稿提出了很多建设性的修改意见和建议，使本研究更加完善，在此一并表示感谢.