孟秋，胡才博，石耀霖

(中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院 中國科學(xué)院計算地球動力學(xué)重點實驗室，北京 100049)( 2019年12月20日收稿； 2020年2月17日收修改稿)

冰凍圈是指地球上存在的所有冰體所構(gòu)成的空間，包括包含冰的云層、地下冰、積雪、水冰、冰川、冰蓋等瞬時、短期、季節(jié)和永久存在的冰體，最早在1923年由Dobrowolsky提出[1]。人類對冰川和冰蓋在調(diào)節(jié)氣候和地球地貌演化等長期效應(yīng)和泥石流、洪水等短期效應(yīng)方面都做了深入的研究[2]。自工業(yè)革命以來，與人類活動密切相關(guān)的全球氣候變化，正在加速冰川消融。研究現(xiàn)代冰川消融乃至消失對地表變形、地應(yīng)力和主要斷層地震危險性的時空變化，具有理論和現(xiàn)實意義。

中國有48 571個冰川，覆蓋了5.18萬km2的地區(qū)，占世界除南極和格陵蘭島之外冰川總和的7.1%[3]。青藏高原地區(qū)擁有世界上除兩極之外數(shù)量最多的冰川，很多冰川是內(nèi)陸河流的源頭[4]，對于水循環(huán)和物質(zhì)平衡具有重要的意義[5]。本文主要關(guān)注祁連山地區(qū)(36°30′～39°30′N,93°30′～103°00′E)的現(xiàn)代冰川。祁連山位于中國西北內(nèi)陸干旱區(qū)，隸屬青藏高原東北緣，氣候?qū)儆跍貛О敫珊祬^(qū)，受大陸性氣候和青藏高原氣候的影響[6-8]。祁連山地區(qū)由許多NW-SE走向的平行山脈和山谷組成，長度約為800 km，寬度約為300 km，東起烏鞘嶺，西止于當(dāng)金山口，向東南延伸到河西走廊，地形從東北向西南逐漸隆升，最高點為團(tuán)結(jié)峰(5 826 m)[3]。

祁連山地塊是青藏高原向北擴(kuò)張的前緣部位，受到青藏高原與北部阿拉善塊體的擠壓，發(fā)生強(qiáng)烈變形[9]，GPS整體呈NE-SW向壓縮，祁連山北緣斷裂滑動主要運動特征為逆沖兼具右旋走滑，呈現(xiàn)出空間不均勻變形的現(xiàn)象[10]。祁連山主要斷裂帶包括祁連山斷裂帶和阿爾金斷裂帶，其中祁連山斷裂帶由多條活動斷裂帶組成，總體為北西西走向伸展，地震活動頻繁[11]。該區(qū)域主要受到擠壓作用，同時也會受邊界走滑作用的影響，形成左旋逆沖兼走滑的構(gòu)造特征[12]。祁連山地震帶被深大走滑活動斷裂包圍，區(qū)內(nèi)發(fā)生過多次大地震，包括1920年海原8.5級地震、1927年古浪8級地震等[13]。

祁連山冰川屬于阿爾卑斯冰川，是冰凍圈的重要組成部分，同時也是內(nèi)陸地區(qū)重要的水源，被稱為“阿爾卑斯固體水庫”，其融水對于中國西北部脆弱的山區(qū)生態(tài)系統(tǒng)起著至關(guān)重要的作用，是中國西北部干旱區(qū)的重要水源之一。同時，也加強(qiáng)了中國地區(qū)與“絲綢之路”其他國家地區(qū)的往來，為中國內(nèi)陸干旱地區(qū)經(jīng)濟(jì)的發(fā)展、文化的繁榮提供了保障[14]。祁連山冰川中，小于1 km2覆蓋面積的冰川數(shù)量最多，1～5 km2的冰川覆蓋面積最大，其中58%面積的冰川位于海拔4 800～5 200 m之間的區(qū)域[3]。

隨著全球氣溫的升高，祁連山地區(qū)的冰川逐漸消融，雪線不斷向后退縮，在1956—2000年間可以達(dá)到419 m/a。1956—2010年期間，祁連山冰川總面積減少420.81 km2，約占總面積的20.88%；體積減少21.63 km3，占總體積的20.26%。其中小型冰川的減少是主要原因，目前海拔4 000 m以下冰川已經(jīng)完全消失，4 500 m以下的冰川跟50 a前相比，消失了一半以上。4 350～5 100 m之間的冰川減小面積占總面積的84.21%。整體看來，冰川消融呈現(xiàn)出東部冰川減少最快、西北部減少較慢的趨勢[3]。

祁連山地區(qū)如此巨量的冰川載荷變化可能會引起地殼變形、地應(yīng)力的時空變化，甚至可能誘發(fā)新的地震活動。歷史上有很多人類活動或者自然環(huán)境變化引起地殼變形、影響地震活動性的先例。Liu等[15]指出，颶風(fēng)可能引發(fā)活動斷裂上的慢地震的發(fā)生；Bettinelli指出喜馬拉雅GPS速度的垂向分量具有季節(jié)性，與恒河平原地表水載荷有關(guān)[16]，并且在冬季小震活動更加頻繁[17]。前人對地下水開采、水庫蓄放水、地?zé)衢_發(fā)高壓注水、油氣注水開采等也進(jìn)行了大量的研究工作。龐亞瑾等[18]模擬計算華北地區(qū)地下水的開采對地殼應(yīng)力的影響，結(jié)果表明地下水的開采有減緩大地震孕育的作用；程惠紅等探討紫坪鋪水庫對汶川地震的觸發(fā)作用，結(jié)果表明，水庫造成的影響恰好處于是否觸發(fā)的邊緣[19-20]，周圍小震應(yīng)該與紫坪鋪蓄水有直接關(guān)系[21]；蘆山地震與水庫蓄水的關(guān)聯(lián)度不大[22]，而卡里巴水庫的蓄水則導(dǎo)致了6.1級地震的發(fā)生[23]。Albaric等[24]研究南澳Paralana地區(qū)大規(guī)模水壓致裂過程與地震釋放的關(guān)系。Zang等[25]研究歐洲地區(qū)地?zé)峁こ虒Φ卣鹞ｋU度的影響；Grogoli等[26]和Kim等[27]則指出，2017年11月韓國5.5級地震則受到地?zé)嵯到y(tǒng)2 a持續(xù)高壓注水的影響，人類工業(yè)活動可能是造成該地震的重要因素；廢水的注入和油氣的開采也可能誘發(fā)地震的發(fā)生[28]，美國Okalahoma地區(qū)油氣開采注水可能造成中小地震數(shù)量增多[29]。在這些問題中，數(shù)值模擬是解決問題的重要手段。

黏彈性有限元模型在當(dāng)今地球動力學(xué)模擬中受到越來越多的重視，在各種地球動力學(xué)過程中得到了廣泛應(yīng)用。Freed和Lin[30-32]利用黏彈性有限元模型研究地震同震和震后應(yīng)力變化，認(rèn)為1971 年San FernandoMW6.7 地震可能觸發(fā)了23 a后1994年NorthridgeMW6.7 地震；Wiseman等[33]利用非均勻黏彈性模型研究2004年蘇門答臘地震由于黏彈性松弛引起的震后變形。沈正康等[34]采用黏彈性分層模型研究東昆侖活動斷裂帶大地震之間的黏彈性觸發(fā)作用，認(rèn)為中下地殼的黏性松弛使庫侖應(yīng)力變化隨時間有增強(qiáng)的趨勢。Hu等[35]建立三維黏彈性有限元模型來研究1960年智利大地震的同震及震后變形，模擬結(jié)果較好地解釋了大地震后35 a的GPS觀測結(jié)果。Luo和Liu[36]建立西藏東部及其鄰區(qū)的三維黏彈塑性有限元模型，考慮重力、構(gòu)造加載、下地殼和上地幔應(yīng)力松弛、斷層弱化等因素的影響，并模擬該區(qū)域歷史大地震序列的時空演化，計算2008年汶川大地震的同震變形及同震庫侖應(yīng)力變化，探討2008年汶川大地震的動力學(xué)成因。不同學(xué)者利用黏彈性力學(xué)模型研究了冰川、地表水、地殼變形和地震之間存在復(fù)雜的相互作用，定量分析冰川載荷變化對地表變形和古今地震活動性變化的影響[37-41]。

本文利用自主開發(fā)的Maxwell黏彈性有限元程序，研究1956—2106年間，在冰川不斷消融引起的地表載荷和兩側(cè)地塊的不斷構(gòu)造擠壓作用下，祁連山地區(qū)地表變形、地應(yīng)力的時空演化，定量評估研究區(qū)域主要斷層帶的地震危險性的時空變化。

1 地質(zhì)背景與有限元模型

為了研究由于冰川的消融及近南北向擠壓對祁連山地區(qū)地表變形、地應(yīng)力及主要斷層帶地震危險性變化的影響，我們建立了祁連山地區(qū)一條典型地質(zhì)剖面的二維平面應(yīng)變Maxwell黏彈性有限元模型[42-43]。剖面左端點為(96°30′E,37°N)，右端點為(98°30′E，39°30′N)，對應(yīng)于圖1(a)中紅色實線位置。

圖1 祁連山地質(zhì)背景與剖面地形圖[46]Fig.1 Geological background and topographic section map of Qilian Mountains

根據(jù)地震層析成像結(jié)果(CRUST1.0)建立有限元模型如圖2(a)所示，模型主要由3層介質(zhì)組成，地表是根據(jù)實際地形插值繪制[44]。莫霍面起伏參考地震層析成像的結(jié)果(CRUST1.0)，莫霍面以上2層為各向同性彈性介質(zhì)，莫霍面以下為各向同性Maxwell黏彈性介質(zhì)。為了突出斷層帶介質(zhì)的差異，地質(zhì)剖面上的5條斷層看作寬度為300 m的斷層帶，斷層的位置和幾何產(chǎn)狀來自前人的地質(zhì)考察結(jié)果[44]，視為橫觀各向同性介質(zhì)[45]。斷層帶介質(zhì)比圍巖要小一些，斷層帶介質(zhì)的G2為周圍介質(zhì)的一半。介質(zhì)材料參數(shù)見表1[46]。本文的有限元模型采用四邊形網(wǎng)格剖分，共有36 096個節(jié)點，35 700個單元，網(wǎng)格剖分見圖2(b)。

圖2 有限元模型與網(wǎng)格剖分圖Fig.2 Finite element model and mesh model

表1 黏彈性有限元模型的參數(shù)表[46]Table 1 Material parameters of the viscoelasticfinite element model

邊界條件：通過GPS插值得出剖面兩側(cè)的GPS速度值，將右側(cè)邊界設(shè)為水平方向位移分量為零，垂直方向剪應(yīng)力為零；左側(cè)邊界水平速度為左右兩側(cè)GPS水平速度之差，垂直方向剪應(yīng)力為零；下邊界水平方向剪應(yīng)力為零，垂直方向位移分量為零；上邊界施加隨時間變化的冰川融化卸載的載荷，相當(dāng)于在上邊界施加了向上的面力，其大小隨冰川空間部位和融化時間而變化。值得注意的是，1956年祁連山地區(qū)的冰川范圍比圖1(a)所示的范圍要大很多。根據(jù)相關(guān)文獻(xiàn)，本文假設(shè)冰川載荷的厚度和分布均隨時間不斷變化：冰川最低海拔隨時間按10 m/a的速度上升，冰川厚度在4 500 m處的消融速率為0.32 m/a[47]，從最高點到冰川最低海拔按線性分布，滿足以下公式

(1)

其中：H為當(dāng)前海拔冰川厚度，H0為海拔4 500 m(h0)處冰川厚度，h為當(dāng)前海拔高度，hlmin為冰川最低海拔。

計算時間階段：1956—2106年，共150 a，計算的時間步長為1 a。由于Maxwell黏彈性模型是線性的，本文只考慮邊界構(gòu)造加載和地表冰川載荷變化的影響，不考慮初始應(yīng)力場的影響，因此在模型中不考慮重力的因素。

2 模擬計算結(jié)果

通過與1956年的初始位移場相減，得到1956—2106年間共150 a的位移應(yīng)力變化量。圖3給出第50、第100和第150年位移的時空變化。1956—2106年間，隨著時間推移，水平位移不斷增加，主要受到模型兩側(cè)擠壓變形的影響，冰川載荷的變化對其影響相對較小(圖3(a))；冰川載荷的影響主要表現(xiàn)在垂直位移上，在冰川載荷的逐漸消融和兩側(cè)擠壓的共同作用下，地表呈現(xiàn)明顯的向上抬升的現(xiàn)象，最大抬升量達(dá)0.55 m以上，主要位于海拔超過4 500 m的高海拔地區(qū)，海拔越高的地區(qū)，冰川載荷減小越多，抬升幅度越大，模型北段的抬升量不夠明顯(圖3(b))。

圖3 1956—2106年間位移變化等值線圖(等值線單位：m)Fig.3 Displacement contour map during 1956-2106

1956—2106年間，水平應(yīng)力、垂直應(yīng)力和剪應(yīng)力都發(fā)生了顯著變化(圖4)。對于水平應(yīng)力和垂直應(yīng)力的變化，大于零代表拉應(yīng)力增加，壓應(yīng)力減小，對應(yīng)暖色區(qū)(紅色)；小于零代表拉應(yīng)力減小，壓應(yīng)力增加，對應(yīng)冷色區(qū)(藍(lán)色)，以下同。

圖4 1956—2106年間應(yīng)力變化等值線圖(等值線單位：Pa)Fig.4 Horizontal stress contour map during 1956-2106

在構(gòu)造擠壓和冰川卸載的雙重作用下，各主要斷層帶的正應(yīng)力和剪應(yīng)力呈現(xiàn)不同的演化趨勢，最終會對地震危險性的改變帶來不同影響，可以通過庫侖應(yīng)力變化進(jìn)行展示。庫侖應(yīng)力變化的計算公式如下

(2)

圖5給出各主要斷層帶中間深度點處庫侖應(yīng)力變化ΔCFS在1956—2106年間隨時間的變化曲線，可以清楚看出位于地質(zhì)剖面不同部位的不同產(chǎn)狀的各條斷層，其庫侖應(yīng)力變化ΔCFS隨時間的變化曲線差異性很大，也意味著不同斷層帶的地震危險性隨時間的變化的差異性很大。祁連山北段斷裂、哈拉湖南山斷裂和哈拉湖盆地北緣斷裂的庫侖應(yīng)力變化ΔCFS在這150 a期間都是大于零的，并且持續(xù)增加，地震危險性也持續(xù)增加，需要重點關(guān)注。右側(cè)2條斷層、特別是祁連山北斷層ΔCFS基本變化不大，增長緩慢，非常接近零，可能反映僅僅在構(gòu)造應(yīng)力作用下，地震危險性增加的速率相對較小，而冰川的融化，特別在冰川融化量最大的最近50 a前后，對融化地區(qū)下方的斷層地震危險性有較大影響。

圖5 各主要斷層帶中間深度點ΔCFS變化曲線Fig.5 Curves of ΔCFS with time at the center of major faults

3 討論

本文研究祁連山地區(qū)現(xiàn)代冰川融化對地表變形、地應(yīng)力和斷層地震危險性時空變化的影響?，F(xiàn)代冰川的分布如圖1所示，文獻(xiàn)表明其正在加速融化[3]。根據(jù)計算結(jié)果，繪制模擬計算50、100、150 a和現(xiàn)今地表抬升速度圖，如圖6(a)所示。圖像表明，當(dāng)?shù)乇泶嬖诒〞r，海拔越高，冰川融化速度越快，地表抬升速度越大。冰川融化速度越快的地區(qū)，地表抬升速度越大；當(dāng)冰川融化后，地表抬升速度依賴于邊界構(gòu)造的擠壓作用。2006年祁連山地區(qū)研究剖面的垂直抬升速率為1～5 mm/a，2009年抬升速率略有減小，這些模擬結(jié)果與前人觀測得到祁連山地區(qū)垂直抬升速率[50-51]較為一致。2056和2106年，冰川完全融化，抬升速率趨于穩(wěn)定，最大值只有1 mm/a左右。

然后對研究區(qū)域的地震震級完備性進(jìn)行驗證。結(jié)果表明，圖1(a)所示的研究區(qū)域內(nèi)1970—2016年間3.0～3.9級地震的數(shù)目為437，4.0～4.9級地震數(shù)目為254，5.0～5.9級地震數(shù)目為67，6.0～6.9級地震數(shù)目為13，7.0～7.9級地震數(shù)目為1，3級以上的地震基本滿足古登堡和里克特地震震級頻度G-R關(guān)系，3.5級以上的地震完備性較好，如圖6(b)所示。圖6(c)給出本文研究的祁連山剖面兩側(cè)30 km以內(nèi)的1956—2016年的3級以上歷史地震分布，這60 a的歷史地震每隔10 a以不同顏色的圓圈表示。從空間上看，在剖面的南段，歷史地震分布最多，主要位于剖面左一和左二的大柴旦—宗務(wù)隆山、哈拉湖南山2條斷層之間，震源深度主要位于5～30 km。研究剖面的中段地震極少，北段有零星的地震分布。這種空間分布是地質(zhì)構(gòu)造條件決定的，不是本文討論的范圍。我們關(guān)心的是從時間來看，1996年以后的地震活動有顯著增強(qiáng)，雖然不能肯定它們就是冰川融化造成的應(yīng)力變化導(dǎo)致的，但這種地震活動性的時間變化與本文計算得到剖面左側(cè)3條斷層地震危險性增加的結(jié)果的確是較為一致的。

圖6 地表抬升速率與地震數(shù)據(jù)分析圖Fig.6 Analysis of surface uplift rates and seismic data

現(xiàn)代冰川的消融，導(dǎo)致本文研究剖面除左側(cè)第1條斷層之外的4條斷層(從南到北依次為：HSF，哈拉湖南山斷裂；HNF，哈拉湖盆地北緣斷裂；YAF，陰凹槽斷裂；QNF，祁連山北緣斷裂)的斷層面正應(yīng)力增加(正應(yīng)力以拉為正)，不僅有利于這4條斷層逆斷層地震的觸發(fā)，而且有利于這4條斷層發(fā)生左旋走滑斷層地震的觸發(fā)。祁連山地區(qū)的逆斷層地震和左旋走滑地震都可能因為現(xiàn)代冰川融化而加速觸發(fā)。

為了解現(xiàn)代冰川融化對祁連山地區(qū)地震危險性變化影響的程度，不妨比較一下歷史大地震對祁連山地區(qū)主要斷裂造成的庫侖應(yīng)力大小。2008年在青藏高原東緣和四川盆地交界的龍門山斷層發(fā)生MW7.9汶川大地震，汶川地震震中離祁連山地區(qū)的距離有800 km以上，造成祁連山地區(qū)主要斷層的庫侖應(yīng)力變化接近于零，不易直接觸發(fā)地震[52]。張貝等[53]利用全球分層介質(zhì)研究2015年4月25日MS8.1尼泊爾大地震的同震應(yīng)力觸發(fā)，對1 000 km以外的祁連山地區(qū)的應(yīng)力變化影響很小。任天翔等[54]利用三維含地形高程的橫向不均勻性橢球型地球有限元模型研究2001年昆侖山口西MS8.1大地震對周圍斷層的同震庫侖應(yīng)力變化的影響，結(jié)果表明祁連山地區(qū)的祁連山斷裂的平均庫侖應(yīng)力變化為-0.001 kPa，宗務(wù)隆—南祁連沖斷裂的平均庫侖應(yīng)力變化為-1.29 kPa，均不利于觸發(fā)地震。通過這些比較，可以認(rèn)識到本文計算的冰川融化與構(gòu)造應(yīng)力疊加對祁連山主要斷層造成的庫侖應(yīng)力變化速率，最大可達(dá)+1 kPa/a以上，而且在150 a內(nèi)不斷累積，這是相當(dāng)可觀的效應(yīng)，不容忽視。

4 結(jié)論

本文利用自主開發(fā)有限元程序，構(gòu)建祁連山地區(qū)的二維平面應(yīng)變黏彈性有限元模型，幾何模型考慮了實際地形和莫霍面起伏，模擬祁連山地區(qū)剖面上邊界構(gòu)造加載和冰川載荷不斷消融對地表變形、地應(yīng)力和斷層地震危險性的時空變化，得出以下初步認(rèn)識：

1)在構(gòu)造加載和冰川持續(xù)卸載的情況下，冰后回彈效應(yīng)明顯，最大垂直抬升可達(dá)0.55 m以上，海拔越高，位移回彈量也越大。隨冰川載荷的減少，地表抬升幅度逐漸增大，海拔較高的中段抬升量最大。

2)1956—2106年間，研究區(qū)域的水平應(yīng)力、垂直應(yīng)力及剪應(yīng)力在冰川載荷卸載和邊界構(gòu)造加載的雙重作用發(fā)生顯著變化；各主要斷層帶的地震危險性呈現(xiàn)差異性的變化趨勢，其中大柴旦—宗務(wù)隆山斷層、哈拉湖南山和哈拉湖盆地北緣斷層ΔCFS大于零，且隨時間一直增大，地震危險性持續(xù)增加，應(yīng)該引起足夠重視；陰凹槽、祁連山北緣2條斷層ΔCFS變化不大，未來發(fā)生地震的風(fēng)險相對較小。

本文只考慮了二維平面應(yīng)變黏彈性模型，未考慮祁連山地區(qū)現(xiàn)代冰川的三維空間展布，對冰川載荷的時空變化估計還比較粗略，下一步將開展祁連山地區(qū)三維黏彈性有限元模型的工作，以期進(jìn)一步研究冰川消融對該地區(qū)地震危險性時空變化的影響。

感謝中國地震局地質(zhì)研究所張竹琪博士在地震震源數(shù)據(jù)和繪圖方面提供的幫助。