東昆侖斷裂的震后粘彈性松弛效應對地殼變形的影響

付譽超 刁法啟 朱亞戈 陳 飛 熊 熊

1 中國地質大學(武漢)地球物理與空間信息學院湖北省地球內部多尺度成像重點實驗室，武漢市魯磨路388號，430074

近30年來，以GPS和InSAR為代表的空間大地測量技術發展迅速，實現了高精度地殼變形監測，可為深入了解不同時空尺度下的地殼變形機理提供基礎資料[1]。某一時段內觀測到的地殼變形是區域地球動力學過程的綜合反映，不僅包含穩態的構造變形，還可能包含強震引起的瞬態同震變形和幾年至幾十年時間尺度的震后變形等其他局部因素[2]。Calais等[3]研究發現，蒙古地區2個8級左右地震在震后約100 a仍能造成每年數mm的地殼變形。震后變形是震后過程研究中不可或缺的約束信息，但對穩態構造運動研究(如斷層穩態滑動速率和塊體運動學特征)而言，震后變形則被視為“噪聲”，需在地殼變形研究中予以考慮[4]。

圖1 研究區構造背景Fig.1 Tectonic settings of the study area

東昆侖斷裂是橫跨青藏高原的大型左旋走滑斷裂之一，沿東西走向展布近2 000 km，以約10 mm/a的運動速率調節南側巴顏喀拉塊體和北側柴達木盆地的差異[5]。自20世紀初以來，東昆侖斷裂帶附近共發生6次7級以上地震(圖1)，分別為1937年花石峽7.5級地震、1963年都蘭7.1級地震、1973年瑪尼7.3級地震、1997年瑪尼7.5級地震、2001年昆侖山口西8.1級地震和2021年瑪多7.4級地震[6-10]。其中，發生時間較近的1997年瑪尼7.5級地震、2001年昆侖山口西8.1級地震和2021年瑪多7.4級地震均觀測到明顯的震后變形[11-14]。有學者分析震后變形的力學機制發現，粘彈性松弛效應是影響長時間尺度下大范圍變形的主要機制，而震后余滑主要影響短期、近場震后變形[15]。上述研究有助于深入了解青藏高原東北緣下地殼/上地幔的流變結構和震后變形過程。然而，前人通常僅針對某一時期的單個震例進行研究，缺乏整條斷層長時間尺度震后效應的綜合分析。

基于此，本文以東昆侖斷裂為例，結合前人對區域流變結構的認識，采用較為合理的介質模型和震源參數對該區域近百年來發生的6個強震的震后松弛效應進行數值模擬。在此基礎上給出該效應對區域地殼變形的影響，并對結果進行分析和討論，從而為其他大型斷裂上強震的震后變形綜合研究提供理論依據。

1 模型構建和參數設置

本文采用分層介質模型，結合地質、歷史地震記錄等資料，考慮較為真實的地震破裂分布情況，基于前人研究得到的區域下地殼/上地幔粘滯系數，綜合分析東昆侖斷裂帶附近區域近百年來發生的6個M≥7地震的震后粘彈性松弛效應造成的地殼變形情況及時空演化過程。

本文采用PSGRN/PSCMP軟件包[16]模擬東昆侖斷裂附近區域強震的震后粘彈性松弛效應。分層介質模型包含彈性的上地殼、粘彈性的下地殼和上地幔。介質模型的彈性參數取自區域地震學研究結果(圖2(a)、(b))，其中下地殼和上地幔的粘彈介質采用Burgers體流變結構表示(圖3)。Burgers體由一個Maxwell體和一個Kelvin體串聯而成，可描述介質的不同松弛時間特征，其本構關系可表示為[17]：

(1)

式中，τ1=η1/μ1，τ2=η2/μ2，τ12=η1/μ2。

在常應力加載狀態下，應變隨時間的變化可表示為：

(2)

式中，應變為瞬態的彈性響應、呈指數衰減的短期粘彈響應和線性變化的長期穩態粘彈響應等3個響應的線性疊加。

下地殼/上地幔的粘滯系數較難確定，但該區域的震后變形研究可為該參數的確定提供約束(圖2(c))。Diao等[18]根據1999～2015年的GPS觀測數據得出該地區下地殼的長期穩態粘滯系數約為1×1018Pa·s。Diao等[15]考慮震后余滑和震后粘彈性松弛效應的綜合影響后，根據2001年昆侖山口西8.1級地震震后4個月的觀測資料發現，下地殼和上地幔的粘滯系數約為1×1018Pa·s。Ryder等[19]以昆侖山口西8.1級地震震后1 a的GPS資料和震后2～5 a的InSAR資料為約束，得到其下地殼瞬態和穩態粘滯系數分別為1×1018Pa·s和1×1019Pa·s。根據上述研究結果，本文設置該地區下地殼瞬態和穩態粘滯系數的下邊界分別為1×1017Pa·s和1×1018Pa·s，對應的上邊界分別為2.5×1018Pa·s和2.5×1019Pa·s(圖2(c))。根據文獻[13]的研究結果，本文將該地區上地幔瞬態和穩態粘滯系數分別設置為1×1018Pa·s和1×1019Pa·s。

VP和VS分別為地震P波和S波波速，ρ為地球介質密度，不同數字序號為不同研究得到的結果，本研究中下地殼穩態粘滯系數η2取1×1018～2.5×1019 Pa·s(淺藍色陰影部分)，橙色五角星為該范圍的平均值圖2 介質模型參數Fig.2 Parameters of the media model

η1和η2分別為瞬態和穩態粘滯系數，μ1和μ2分別為2個彈簧的剪切模量圖3 Burgers體結構Fig.3 Burgers body

本文根據2001年昆侖山口西8.1級地震和1997年瑪尼7.5級地震的同震破裂分布深度[20]以及該地區中小地震分布深度(圖2(d))，綜合確定彈性上地殼的厚度[21]。大地測量數據反演結果表明，上述2個地震的同震破裂主要分布在20 km以上，且95%區域的中小地震也分布在20 km以上。上述統計結果與Ryder等[19]選取的該地區彈性上地殼的厚度(16.5 km)較為接近，因此本文將彈性上地殼厚度設置為16.5 km。本文進一步根據區域地震學成像結果[22]和震后變形研究結果[19]，將下地殼厚度設置為65 km。

此外，同震破裂模型也會影響粘彈性松弛效應。本文采用由InSAR資料反演的精細位錯模型[14,20,23]對1997年瑪尼7.5級地震、2001年昆侖山口西8.1級地震和2021年瑪多7.4級地震進行研究，而對早期發生的3個地震采用Xiong等[24]估算的結果進行研究。

2 結果與分析

由于該地區下地殼粘滯系數存在較大的不確定性，會對結果產生顯著影響，因此本文將粘滯系數設定為圖2(c)中取值范圍的平均值(穩態粘滯系數設置為1.3×1019Pa·s)進行模擬分析。由圖4(a)可知，由于1937年花石峽7.5級地震、1963年都蘭7.1級地震和1973年瑪尼7.3級地震的離逝時間相對較長[24]，因此這些地震的震后粘彈性松弛造成的地表變形相對較弱(<1 mm/a)。說明對于現今地表而言，該區域的震后變形主要由離逝時間較短的2001年昆侖山口西8.1級地震和1997年瑪尼7.5級地震[24]共同構成。模擬結果顯示，2010～2020年震后粘彈性松弛效應使得東昆侖斷裂帶兩側的地表顯著變形，其中最大地表變形發生在昆侖山口西8.1級地震震源區南側，其東向變形速率可達7.6 mm/a，而斷裂帶北側的西向最大地表變形速率約為6.0 mm/a。此外，震后粘彈效應在1997年瑪尼7.5級地震震源區附近引起的最大地表變形速率約為6.0 mm/a。

對青?，敹?.4級地震造成的震后粘彈效應進行模擬，由圖4(b)可見，2021～2030年該地區預計可產生最大約3.9 mm/a的地表變形。該結果可為相關震后變形研究提供理論依據，同時表明在進行其他研究時，應考慮地表變形的影響。

為更直觀地展示東昆侖斷裂附近強震震后松弛效應的影響，將模擬的震后位移場投影到平行于斷層走向的剖面上，并將其與震間變形模型預測的地表變形速率進行比較。圖4(c)、(d)分別將速度場投影到A-A′和B-B′剖面上，其中剖面A-A′ 和B-B′處的斷層滑動速率(10 mm/a)取自地質學和大地測量學研究結果[4,18]，斷層閉鎖深度(16.5 km)則取自Ryder等[19]的研究結果。

由圖4可見，剖面A-A′上變形速率分量主要反映的是1997年瑪尼7.5級地震和2001年昆侖山口西8.1級地震的震后變形結果。由圖4(c)可見，2005～2010年和2010～2020年震后粘彈性松弛效應對剖面A-A′造成的最大影響分別約為27.2 mm/a和8.4 mm/a，均高于震間變形模型預測的變形速率。對瑪多7.4級地震而言，2021～2025年和2025～2030年粘彈性松弛效應對剖面B-B′造成的最大影響分別約為3.4 mm/a和2.6 mm/a，其影響不可忽略。

模擬結果顯示，距斷層較近區域(≤100 km)的震后松弛效應的影響與震間變形模型預測的結果大致相同(圖4(c)、4(d))，說明若僅根據觀測到的地表速度場來研究斷層的滑移率而不考慮其震后松弛效應，必然會使研究結果產生較大偏差。

3 討 論

3.1 粘滯性參數對研究結果的影響

3.1.1 下地殼粘滯系數的影響

研究表明[26]，震后粘彈性松弛效應響應時間主要受粘滯系數約束。本文以2001年昆侖山口西8.1級地震為例，選取震源區附近一點(92.53°E, 36.33°N)，討論不同下地殼粘滯系數的震后地表水平位移場時空變化特征(圖5)。由圖可見，下地殼粘滯系數變化對粘彈性松弛效應的弛豫時間影響顯著。粘滯系數越小，早期地表變形速率越快，粘彈性松弛效應衰減也越快；粘滯系數越大，應力松弛越慢，地表變形速率越慢，震后粘彈性松弛效應衰減越慢，持續時間也更長。

基于此，為研究粘滯系數不確定性對模擬結果的影響，對前文選取的下地殼粘滯系數范圍進行討論(圖2(c))。由圖6(a)可見，2005～2010年震后松弛效應在斷層附近區域造成的地殼變形速率為9.6～44.4 mm/a；由圖6(b)可見，2010～2020年變形速率迅速衰減至mm量級，但仍可達3.6～9.0 mm/a；由圖6(c)可見，2021～2025年東昆侖斷裂附近系列強震的震后松弛效應造成的地表變形速率可達6.2～11.7 mm/a；由圖6(d)可見，2025～2030年該效應引起的周邊地表變形速率預計約為5.9～7.6 mm/a。上述分析表明，當下地殼粘滯系數在取值范圍內變化時，以1997年瑪尼7.5級地震、2001年昆侖山口西8.1級地震和2021年瑪多7.4級地震等為代表的東昆侖斷裂帶附近強震的震后松弛效應，在震后數10 a內對震源區附近地表造成的影響仍可達每年數mm。

(a)和(b)表示前5次地震的影響；(c)和(d)表示所有6次地震的影響圖6 下地殼粘滯系數對地表變形的影響Fig.6 Effect of lower crust viscosity on surface deformation

由圖7(a)可見，2005～2010年前5個地震的震后粘彈性松弛效應對剖面A-A′造成的影響可達7.0 ～43.0 mm/a，即便是距斷層超過100 km的區域，斷層兩側的變形速率仍大于震間變形模型的模擬結果；由圖7(b)可見，2010～2020年斷層附近A-A′剖面變形速率可達3.4～9.4 mm/a；由圖7(c)可見，對于發生時間較近的2021年瑪多7.4級地震，2021～2025年震后松弛效應對剖面B-B′造成的影響約為2.3～9.5 mm/a；由圖7(d)可見，2025～2030年震后松馳效應對剖面B-B′ 造成的影響預計會減小至1.7～6.2 mm/a，其影響不可忽略。

(a)和(b)表示前5次地震的影響；(c)和(d)表示瑪多7.4級地震的影響；灰色陰影區域表示模擬得到的跨斷層變形速率范圍，黑色曲線和紅色曲線分別表示模擬結果的平均值和根據震間變形模型得到的結果[25]圖7 下地殼粘滯系數對跨斷層滑動速率的影響Fig.7 Effect of lower crust viscosity on slip rate across the faults

3.1.2 下地殼介質非均質性、流變模型及上地幔粘滯系數選取

研究發現，東昆侖斷裂帶南北兩側下地殼粘滯系數存在明顯差異[21,27]，說明建立更為精細的三維流變結構有限元模型可進一步深入理解震后變形過程。然而，目前的觀測數據對該區域下地殼粘滯系數的約束能力較差，不同研究給出的該地區下地殼粘滯系數結果之間的差異(10～100倍)甚至遠高于斷裂帶兩側粘滯系數之間的差異(約10倍)[28]。因此，本文采用橫向均勻的介質模型進行震后粘彈性松弛模擬研究。

相比于Maxwell體，Burgers體需要考慮更多參數，并且能同時兼顧震后瞬態和長期的變形模擬[29-30]。本文在采用Burgers體進行模擬時，根據前人研究結果，將瞬態粘滯系數設定為穩態粘滯系數的1/10[13]。模擬結果表明，由于該區域的震后松弛主要由下地殼的粘性流動引起[19]，上地幔的貢獻相對較小，且現有數據對該區域上地幔粘滯系數的約束能力較弱，研究相對較少，因此本文將上地幔粘滯系數設為固定值。

3.2 震后粘彈性松弛效應對斷層滑動速率估算的影響

Bell等[31]分析發現，1997年瑪尼7.5級地震前瑪尼斷裂周邊的地表變形速率約為1 mm/a，在瑪尼7.5級地震發生10 a后，瑪尼斷裂附近的地殼變形速率可達8 mm/a，遠大于震前的變形速率。由圖4(a)可見，2005～2010年震后松弛效應在瑪尼斷裂附近造成的地表變形速率最大約為6.0 mm/a，說明震后松弛效應或許可以解釋1997年瑪尼7.5級地震附近現今地殼變形速率明顯大于震前變形速率這一現象。Li等[32]研究發現，2009～2015年昆侖山口西8.1級地震震源區域震后地表變形速率相較于震前增加約40～80 mm/a；Liu等[33]研究發現，2007～2015年昆侖山口西8.1級地震的震后變形速率約為10 mm/a。由圖6(a)可見，昆侖山口西8.1級地震震后松弛效應在震后10 a內造成的地表平均變形速率約為9.6～44.4 mm/a，與前人研究結果基本一致。

Gan等[34]利用跨東昆侖斷裂GPS觀測數據和彈性半空間模型得到東昆侖斷裂中段的滑動速率約為20 mm/a，而地質學研究結果[35]與其他大地測量研究結果[5]中的滑動速率分別約為8～12 mm/a和10～15 mm/a。Wang等[36]與Li等[32]將地質學和大地測量研究得到的滑移速率結果差異歸因于不同的空間位置，Gan等[34]推斷該差異可能是由震后粘彈性松弛效應所致。

基于此，本文根據Zheng等[37]發表的GPS速度場，采用2001年昆侖山口西8.1級地震近場GPS臺站數據(圖1中紅色三角形)構建垂直東昆侖斷裂中段的速度剖面，并使用一維模型搜索該段斷層的滑動速率。已有研究表明，垂直于斷層走向的高密度GPS觀測數據可為斷層滑移率估算提供可靠約束[38]。本文采用前人基于彈性位錯理論提出的震間變形模型進行反演[25]，該模型中平行于斷層走向的臺站變形速率V(x)可用斷層閉鎖深度D、滑動速率V和臺站距斷層的距離x來表示：

(3)

圖8 斷層滑移速率反演結果對比Fig.8 Comparison of inverted results of fault slip rates

由式(3)可得到一條關于x的反正切函數曲線，即震間變形曲線(圖8(c)中綠色和黑色虛線)。需要說明的是，本文采用的是Zheng等[37]發表的速度場數據，但并未考慮昆侖山口西8.1級地震震后粘彈性松弛效應的影響，因此分2種情況展開研究：1)情況1，直接采用Zheng等[37]發表的速度場進行斷層滑移速率反演；2)情況2，首先在Zheng等[37]發表的速度場基礎上扣除該時段內模擬的震后粘彈性松弛效應，以此獲得震間穩態的地表速度場，然后進行斷層滑動速率反演。初步測試發現，地表變形數據對斷層閉鎖深度D并不敏感，因此對斷層閉鎖深度D取固定值。由于該斷層閉鎖深度D與區域彈性上地殼層厚度[19]具有較好的一致性，因此在參考昆侖山口西8.1級地震同震破裂深度分布結果[20]時，將斷層閉鎖深度D固定為16.5 km。各臺站原始GPS觀測數據和扣除震后粘彈性松弛效應后的地表變形數據如圖8(c)中藍色和紅色誤差棒所示。由圖8(a)、(b)可見，情況1得到的東昆侖斷裂中段的滑動速率V約為13 mm/a，情況2得到的滑動速率V約為10 mm/a。需要說明的是，由于數據限制，本文未考慮東昆侖斷裂南側塊體內部次級斷層的作用，但這些次級斷層的滑動可能也會對反演結果產生一定的影響[38]。

由圖8可見，情況2反演得到的東昆侖斷裂中段的滑移速率結果與地質學研究結果[35]具有更好的一致性，而情況1反演的結果則會高估約30%的斷層滑移速率。Li等[32]利用2001年昆侖山口西8.1級地震震后GPS觀測數據(2014～2017年)反演得到的東昆侖斷裂中段的滑移速率(17.4±1.2 mm/a)遠高于使用震前GPS數據(1991～2001年)反演結果(7.2±2.2 mm/a)；而Wang等[36]采用1999～2014年GPS數據，剔除該時段內強震的同震和震后變形影響后得到的該段斷層滑移速率(10.0±0.3 mm/a)與本文結果十分吻合，說明在穩態斷層滑移速率研究中，需要考慮并扣除強震震后粘彈性松弛效應。

由圖4(c)可見，對于現今地表而言，震后松弛效應對昆侖山口西8.1級地震破裂帶附近地殼變形的最大影響可達8.4 mm/a，與震間積累的變形速率影響相當。因此，在利用現今大地測量觀測資料研究斷層滑動速率時，震后松弛效應的長期影響不容忽視。

3.3 不同震后變形機制的影響

震后變形過程可能是震后斷層余滑、下地殼/上地幔的粘彈性松弛和震后孔隙流體調整的共同作用結果[39]。Liu等[33]認為，2001年昆侖山口西8.1級地震震后短期變形機制應該包括余滑和粘彈性松弛效應；賀鵬超等[27]認為，忽略震后余滑的影響會低估介質的粘滯系數；Zhao等[28]研究表明，余滑效應會影響下地殼瞬態粘滯系數。由于震后斷層余滑和孔隙流體調整可能僅作用在震后很短時間內，且影響范圍較小，因此在考慮震后長期影響時，需要將震后粘彈性松弛效應視為主要因素。

4 結 語

本文以東昆侖斷裂為例，構建粘彈變形模型，并設定合理的粘滯性結構和地震位錯模型，模擬近百年來該斷裂附近6個M≥7地震的震后粘彈性松弛效應，并分析該效應對穩態地殼變形的影響。結果表明：

1)東昆侖斷裂附近歷史強震的震后粘彈性松弛效應在2010～2020年對地殼造成的最大變形速率約為7.6 mm/a，而該效應造成的最大跨斷層變形速率約為8.4 mm/a，大于震間模型估計的斷層兩側地表變形速率。

2)2021年瑪多7.4級地震震后粘彈效應預計將在2025～2030年造成最大跨斷層變形速率約為2.6 mm/a，若忽略震后粘彈性松弛效應，將會高估震間斷層滑動速率。

3)考慮震后粘彈性松弛效應得到的東昆侖斷裂中段滑動速率與地質學研究結果具有更好的一致性，而不考慮震后粘彈性松弛效應則會高估30%的斷層滑移率。

致謝:2021年瑪多7.4級地震震源機制解和區域歷史地震數據分別取自美國地質調查局USGS(https:∥earthquake.usgs.gov/earthquakes/map)和中國地震臺網中心國家地震科學數據中心(http:∥data.earthquake.cn/index.html)；文中大部分圖件由GMT軟件繪制，在此一并表示感謝。