2021年瑪多MS7.4地震的震后效應(yīng)模擬

談洪波 王嘉沛 楊光亮 陳正松吳桂桔 申重陽 黃金水

1)中國科學(xué)技術(shù)大學(xué)，地球和空間科學(xué)學(xué)院，合肥 230026 2)中國地震局地震研究所，武漢 430071

0 引言

地震發(fā)生后，中國地震局迅速組織實(shí)施了地震科學(xué)考察，開展了多手段聯(lián)合觀測(包括重力和GNSS等)，為震后地震趨勢判定、 地震應(yīng)急救援等工作提供資料。 震源理論的發(fā)展及觀測數(shù)據(jù)的不斷豐富，使我們可在地震發(fā)生后快速確定震源的基本參數(shù)，也使得將地震位錯理論應(yīng)用于震后趨勢分析成為可能。 Zhang等(2014)發(fā)展了一套反演地震破裂過程的迭代反褶積方法，并聯(lián)合高頻GPS和強(qiáng)震數(shù)據(jù)對日本MW9.0地震的破裂過程進(jìn)行了反演。 劉剛等(2015)利用遠(yuǎn)場地震波和GPS數(shù)據(jù)給出了2015年尼泊爾MS8.1地震的破裂過程。 Feng等(2017)利用InSAR和GPS數(shù)據(jù)反演了2015年尼泊爾MW7.8地震及其MW7.2余震的斷層滑動分布。 談洪波等(2013)模擬了玉樹MS7.1地震的同震重力變化并利用絕對重力觀測結(jié)果進(jìn)行了檢驗(yàn)。 Zhang等(2016)對日本MW9.0地震的同震模擬結(jié)果及觀測形變、 重力變化進(jìn)行了比較，發(fā)現(xiàn)兩者符合較好。 Riguzzi等(2019)模擬了均勻彈性半空間中2009—2016年意大利中部4次MW6地震的同震重力效應(yīng)，并就模擬結(jié)果進(jìn)行了合理解釋。 但實(shí)際地球要復(fù)雜得多，地球曲率、 層狀構(gòu)造、 黏滯性、 自重效應(yīng)、 三維構(gòu)造等均會對計算結(jié)果產(chǎn)生不同程度的影響(Pollitz，1996； Sunetal.，2002； 孫文科，2008，談洪波等，2010)。 同震只是介質(zhì)在完全彈性情況下的瞬時效應(yīng)，對于長時間尺度的震后效應(yīng)而言，其形變和重力變化有可能完全改變同震效應(yīng)的形態(tài)，利用震后數(shù)據(jù)獲取斷層和地球內(nèi)部信息也是非常有意義的研究內(nèi)容(張晁軍等，2008； 賀鵬超等，2018； 趙斌等，2020)。 因此，本文將采用更為接近實(shí)際情況的彈性—黏彈分層介質(zhì)模型，并考慮自重的影響，從形變和重力的角度對該地震的同震及其震后效應(yīng)進(jìn)行理論模擬，為利用實(shí)際形變和重力資料解釋此次地震孕震過程(包括同震和震后)的相關(guān)研究提供理論依據(jù)。

1 理論基礎(chǔ)

均勻彈性半無限空間地球模型的位錯理論由Okada(1985)和Okubo(1992)完善和發(fā)展。 假設(shè)地震破裂發(fā)生在一個平面上，斷層如圖 2 所示，長寬分別為L和W，傾角為δ； 斷層底部深度為d； 斷層面上的走滑、 傾滑和張裂方向錯動量分別為U1、U2、U3，在地表(x1，x2，0)處引起的形變Δu和重力變化Δg可分別表示為(Okada，1985； Okubo，1992)

圖1 研究區(qū)主要構(gòu)造分布Fig. 1 Distribution of the main structures in the study area.黑點(diǎn)為截至2021年5月30日的ML>2余震

圖 2 矩形斷層位錯模型示意圖Fig. 2 Dislocation model of a rectangle fault.

圖 3 SLS流變模型Fig. 3 The SLS rheological model.

(1)

Δg(x1，x2)={ρG[U1Sg(ξ，η)+U2Dg(ξ，η)+U3Tg(ξ，η)]+ΔρGU3Cg(ξ，η)}‖-βΔh(x1，x2)

(2)

其中，G為萬有引力常數(shù)；ρ、 Δρ分別為介質(zhì)密度、 張裂紋內(nèi)密度與介質(zhì)密度之差；S(ξ，η)、D(ξ，η)、T(ξ，η)和Sg(ξ，η)、Dg(ξ，η)、Tg(ξ，η)、Cg(ξ，η)為系數(shù)； 自由空氣重力梯度β=0.309×10-5m/s2，Δh為地表高程變化，Chinnery記號‖表示置換(Okada，1985)。

考慮地球自重影響時，Rundle(1980)首次給出了基于無限小靜態(tài)形變理論動力平衡方程：

·Γ+ρψ+u·g-g(ρu)=0

(3)

(4)

Γ=(λ·u)I+μ[u+(u)T]

(5)

式中，λ為拉梅常數(shù)；μ為切變模量；I為單位張量。 根據(jù)Adams-Williamson條件(Longman，1963)，方程(3)可變換為(Wangetal.，2006)

Γ+ρguz-(ρ

(6)

對于黏彈介質(zhì)，一般基于彈性-黏彈線性對應(yīng)原理(Fung，1965)： 黏彈體在頻率域的應(yīng)力-應(yīng)變關(guān)系與彈性體在時間域的一致，故在頻率域中黏彈問題可采用彈性問題的處理方式計算后再通過逆傅里葉變換回到時間域(Matsuuraetal.，1981)。 Wang等(2006)采用SLS(Standard Linear Solid)流變體將彈性體和黏彈性體用統(tǒng)一的模型表示(圖 2)，這時復(fù)切變模量和復(fù)拉梅常數(shù)可表示為

(7)

(8)

聯(lián)合方程(5)—(8)，根據(jù)正則化Haskell(1953)傳播算法，可將邊值問題轉(zhuǎn)變成代數(shù)問題，根據(jù)給定的地球參數(shù)求出頻率域中對應(yīng)的格林函數(shù)(即單位點(diǎn)脈沖位錯響應(yīng))，再采用抗假頻技術(shù)(Kindetal.，1982)對頻率域中的格林函數(shù)進(jìn)行傅里葉逆變換回到時間域，最后通過式(1)和(2)的線性疊加求出地表形變和重力變化( Wang，1999； Wangetal.，2006)。

2 介質(zhì)模型和發(fā)震斷層模型

前期的研究結(jié)果表明地殼分層對地表同震效應(yīng)的影響可達(dá)20%以上(Sunetal.，2002； 孫文科，2008； 談洪波等，2010)，因此需要考慮垂向分層結(jié)構(gòu)。 本文參考Crust1.0(4)https:∥igppweb.ucsd.edu/～gabi/crust1.html。地殼模型，將研究區(qū)地球分為上、 中、 下地殼和上地幔，具體參數(shù)見表1。 此外，也考慮了下地殼和上地幔介質(zhì)的黏滯性，黏滯系數(shù)η參考張晁軍等(2008)、 賀鵬超等(2018)和趙斌等(2020)對青藏高原的研究結(jié)果。

表1 地殼—上地幔介質(zhì)分層結(jié)構(gòu)模型Table1 The layered structure of the crust-upper-mantle based on Crust 1.0

圖4 斷層滑動分布(據(jù)USGS結(jié)果繪制)Fig. 4 The slip distribution on the fault plane(according to USGS results).

3 同震效應(yīng)

圖5 地表同震效應(yīng)Fig. 5 Surface co-seismic effects.a 經(jīng)度方向位移； b 緯度方向位移； c 垂直位移； d 重力變化

4 震后效應(yīng)

同時我們也對瑪多、 稱多和達(dá)日縣城處的震后形變和重力變化趨勢進(jìn)行了模擬(表2，圖 7)，這3個點(diǎn)的震中距分別為36km、 172km和156km。 圖7b為經(jīng)度方向的位移，瑪多縣城同震時向W運(yùn)動了126mm，400a后位移量將達(dá)197mm，增加了71mm(同震效應(yīng)的56%)； 同震時稱多和達(dá)日縣城向E分別運(yùn)動了8mm和19mm，400a后位移量將達(dá)39mm和43mm，分別增加了31mm(同震效應(yīng)的387%)和24mm(同震效應(yīng)的126%)； 整體而言，瑪多縣城震后EW向位移量在300a后基本穩(wěn)定不變，而稱多和達(dá)日縣城在50a后基本穩(wěn)定不變。 圖7c為緯度方向的位移，瑪多和稱多縣城同震時分別向N運(yùn)動了151mm和6mm，400a后位移量將達(dá)198mm和18mm，分別增加了47mm(同震效應(yīng)的31%)和12mm(同震效應(yīng)的200%)； 同震時達(dá)日縣城向S運(yùn)動了22mm，400a后位移量將達(dá)67mm，增加了45mm(同震效應(yīng)的204%)； 瑪多縣城震后SN向位移量在100a后基本穩(wěn)定不變，而稱多和達(dá)日縣城的位移量基本穩(wěn)定所需時間分別為20a和100a。

表2 部分縣城位置的震后效應(yīng)Table2 The post-seismic effects on some of the cities

圖 7d 為垂向位移，同震時瑪多縣城下沉3mm，5a后累計下沉量達(dá)5.6mm，隨后地表開始持續(xù)抬升，400a后抬升量將達(dá)46.5mm，400a間地表抬升量達(dá)49.5mm(同震效應(yīng)的16.5倍)； 稱多縣城在同震時地表抬升了3mm，震后5a間有微弱的下沉，隨后持續(xù)抬升，400a后抬升量達(dá)31.7mm，抬升量增加了28.7mm(同震效應(yīng)的9.5倍)； 而達(dá)日在同震時下沉4.3mm，震后5a間有微弱的抬升，隨后持續(xù)下沉，400a后下沉量達(dá)65.9mm，下沉量增加了61.6mm(同震效應(yīng)的14.3倍)； 3個縣城震后的垂直位移量在300a后基本穩(wěn)定不變。 圖7e為重力變化，與垂直位移震后變化趨勢類似并呈現(xiàn)一定的負(fù)相關(guān)性： 同震時瑪多縣城的重力值增加5.0μGal，5a后累計增加到5.4μGal，隨后開始持續(xù)減小，400a后重力值減小為4.4μGal，400a間重力值減小了9.4μGal； 稱多縣城在同震時地表重力值減小了1μGal，震后5a間有微弱的增加(0.08μGal)，隨后持續(xù)減小，400a后重力值減小為6.5μGal，400a間重力值減小了5.5μGal； 而達(dá)日在同震時的重力值增加2.0μGal，震后5a間有微弱的減小(0.19μGal)，隨后持續(xù)增大，400a后重力值增大為13.5μGal，400a間重力值增大了11.5μGal； 3個縣城震后地表重力變化在300a后基本穩(wěn)定不變。

圖7 震后長期時變效應(yīng)Fig. 7 Long-term time-varying effects after the earthquake.a 所選點(diǎn)與斷層的相對位置； b 經(jīng)度方向位移； c 緯度方向位移； d 垂向位移； e 重力變化

5 討論與結(jié)論

文中基于彈性-黏彈分層半空間中的矩形位錯理論，采用USGS發(fā)布的由遠(yuǎn)場地震波反演的瑪多MS7.4地震發(fā)震斷層模型，從理論上模擬了此次地震引起的地表同震及震后形變和重力變化，通過分析可得到如下主要結(jié)論：

圖8 震后400a間的黏彈性松弛效應(yīng)與同震效應(yīng)比值分布Fig. 8 The distribution of the ratios between the viscoelastic relaxation effects during 400 years after the Maduo earthquake and the coseismic effects.a 經(jīng)度方向位移； b 緯度方向位移； c 垂向位移； d 重力變化

圖9 震后5a間的黏彈性松弛效應(yīng)與同震效應(yīng)比值分布Fig. 9 The distribution of the ratios between the viscoelastic relaxation effects during 5 years after the Maduo earthquake and the coseismic effects.a 經(jīng)度方向位移； b 緯度方向位移； c 垂向位移； d 重力變化

圖10 同震模擬結(jié)果與GNSS觀測結(jié)果的比較Fig. 10 Comparison between the coseismic simulation and GNSS observation results.藍(lán)色和紅色箭頭為GNSS觀測結(jié)果，紫色和綠色為模擬結(jié)果

致謝USGS提供了發(fā)震斷層模型，中國地震局地震研究所趙斌博士提供了同震GNSS數(shù)據(jù)，德國波茨坦地球科學(xué)中心汪榮江教授提供了PSGRN/PSCMP程序； 審稿專家為本文提供了寶貴建議。 在此一并表示感謝!