2021年5月22日青?，敹郙S7.4地震同震變形的數值模擬及其動力學啟示*

孟思晨,孟秋,陳啟志,胡才博

(中國科學院大學地球與行星科學學院 中國科學院計算地球動力學重點實驗室, 北京 100049) (2021年11月14日收稿; 2022年4月26日收修改稿)

2021年5月22日2時4分,青海果洛州瑪多縣(34.59°N,98.34°E)發生MS7.4地震(圖1),震源深度17 km(中國地震臺網正式測定)[1]。該地震極震區的烈度最大值為Ⅹ度,Ⅵ度及以上區域面積約為53 704 km2,宏觀震中和微觀震中分別為瑪多縣瑪查理鎮和黃河鄉[2]。從發震構造上看,青?，敹嗟卣鸢l生在青藏高原中部的Ⅱ級地塊巴顏喀拉塊體北部,距其北邊界東昆侖斷裂帶約70 km。近年來,巴顏喀拉塊體周緣發生一系列強震,如2008年汶川8.0級地震、2011年玉樹7.1級地震、2013年蘆山7.0級地震、2017年九寨溝7.0級地震,而瑪多地震卻是唯一發生在巴顏喀拉塊體內部的強震[3],因此受到科研人員的廣泛關注。人們從震源機制、震源破裂過程、發震斷裂、余震序列、同震滑動分布、震區構造應力場、庫侖應力變化和地震危險性等方面對瑪多地震進行了研究。

F1：達日斷裂帶;F2：瑪多—甘德斷裂帶;F3：東昆侖斷裂帶;F4：昆中斷裂帶;F5：鄂拉山斷裂帶,斷層數據引自詹艷等[4]。圖1 瑪多地震的構造背景和地震活動性Fig.1 The tectonic setting of the Maduo earthquake and the historical seismicity in research area

中國地震臺網、美國地質調查局等機構對此次地震的震源機制和震源破裂過程反演后認為該地震為一次左旋走滑型地震事件,兼有正斷分量[4],為雙側破裂方式,震區構造應力場最大主應力軸方向為NE-SW,最小主應力軸方向為NNW-SSE[1]。李志才等[5]對GNSS數據進行反演確定的地震矩、斷裂的幾何參數及同震滑動分布,也揭示了發震斷裂左旋走滑特性,同時得到同震滑動量的數值及分布。余震雙差重定位研究發現,余震帶與昆侖山口—江錯斷裂距離最近且兩者走向有一定的重合度,說明瑪多地震的發震斷裂為昆侖山口—江錯斷裂[3,6],與野外地震考察結果相一致[7-10]。前人對于瑪多地震同震及震后變形研究可以分為大地測量數據反演及同震、震后變形及相關應力變化模擬多個方面。研究者基于InSAR、GPS等大地測量數據與圖像,揭示出最大滑動量為6.0 m[11-12],并確定了同震變形、靜態滑動分布以及斷層破裂過程[13],量化了地震間應變率、斷層幾何以及同震滑動分布[14]?；贗nSAR數據反演的斷層滑動模型和PREM地球模型,研究者計算了此次地震的理論同震位移、大地水準面、重力及應變的變化,得到地殼變形的空間分布,為研究地震斷層及解釋大地測量數據提供了參考[15]。在同震、震后變形模擬及庫侖應力計算方面,楊光遠等[16]使用Coulomb3.3程序[17-18]計算瑪多地震引起的地震破裂面及周邊主要斷層上的庫侖應力變化;談洪波等[19]和Li等[20]使用均勻半空間解析解PSGRN/PSCRM程序[21]對瑪多地震同震、震后變形、重力變化及發震斷層及周圍區域的庫侖應力變化進行了計算模擬。

前人使用解析法或半解析法對同震、震后及庫侖應力變化進行理論計算,模型比較簡單,計算速度快,因此,在大地震發生不久,通常會有相關的同震、震后及庫侖應力變化成果發布,對大地震周圍區域的地震危險性做出第一時間的評估和判斷。發震斷層周圍介質的非均勻性、幾何的復雜性,對大地震的同震變形影響很大,簡單的解析解和半解析解很難考慮這些復雜因素,而數值模擬可以發揮獨特的優勢。在數值模擬中,有限元方法可以根據實際地質體的幾何、材料及邊界條件靈活地進行網格劃分,同時進行并行計算,實現對復雜地質體高分辨率模擬。分裂節點法通過將位錯加入到載荷向量而不改變剛度矩陣,可以簡潔有效地將斷層位錯引入有限元計算。為此,本文采用C語言編寫了分裂節點法的三維并行彈性有限元程序,用該程序計算了青?，敹嗟卣鸬耐鹱冃魏蛻ψ兓?并對研究區域主要斷層地震危險性進行探討。

1 有限元方法

1.1 使用分裂節點法的有限元模型構建

用有限元方法模擬地震變形時,首先面臨的問題是如何將斷層的同震位錯引入有限元數值計算中。分裂節點法可以簡單有效地解決該問題[22]。如圖2所示,單元1、2 分別位于斷層的左右兩側,U為位移,上標表示單元號,下標表示單元內的節點號;ΔU表示位錯量[22]。

圖2 分裂節點法的示意圖Fig.2 Schematic diagram of split node technique

則有

(1)

(2)

(3)

(4)

單元1的單元剛度矩陣方程為

(5)

同樣,單元2的單元剛度矩陣方程為

(6)

總體剛度矩陣方程為

(7)

由此可以看出,在分裂節點法中,斷層位錯的引入體現在單元、總體的載荷向量上,而對于單元、總體的剛度矩陣是沒有影響的,大大簡化了剛度矩陣和載荷向量的組裝,且對于計算的準確性和穩定性沒有影響[19],體現了該方法的簡潔性和有效性。

1.2 有限元與理論解結果對比

為驗證我們自主研發的采用分裂節點技術的三維并行彈性有限元程序的有效性,我們建立了一個三維走滑斷層地震模型(模型1)和一個三維逆斷層地震模型(模型2),分別利用有限元程序和基于半無限空間位錯理論[23]的Coulomb 3.3[17-18]及edgrn/edcmp[21]計算2個斷層模型的地表同震位移場,并對比兩者的結果。2個斷層模型的參數如表1所示。

表1 走滑和逆沖斷層地震模型的材料參數與斷層參數Table 1 The material parameters and fault parameters of strike-slip and reverse fault earthquakes

圖3分別展示了理論解和有限元計算得到的走滑斷層和逆斷層的同震地表位移的等值線圖和剖面結果圖。圖3(a)～3(f)左側為Coulomb3.3[17-18]解析解結果,右側為有限元模擬結果,表明有限元得到的走滑斷層地震和逆斷層地震引起的地表同震位移的分布形態與理論解有很好的一致性。圖3(g)～3(i)分別給出了有限元和理論解的過斷層中心點垂直斷層走向的剖面上走滑斷層地震的南北向同震位移、逆斷層地震的東西向位移和垂直向位移的對比圖,FEM表示有限元模擬結果,Theory1表示Coulomb3.3[17-18]計算結果,Theory2表示edgrn/edcmp[21]計算結果,表明有限元與理論解有很好的一致性。

圖3 走滑斷層和逆斷層地震引起的同震位移的有限元模擬結果與解析解對比Fig.3 Comparison of the analytical solutions with finite element models simulation results for the coseismic displacement components of the strike-slip fault and thrust fault

2 青?，敹嗟卣鸬耐鹱冃魏蛻ψ兓瘮抵的M與比較

2.1 有限元模型的設置

采用前人反演得到的非均勻位錯模型以及我們基于分裂節點技術[22]的有限元方法,本文計算了瑪多地震的同震位移和同震應力場。有限元模型東西長1 100 km,南北長1 000 km,深度100 km,共有64個并行計算分區,有限元網格數量為2 171 715,可以實現較高分辨率的有限元并行計算。本文采用中國地震局地質研究所的非均勻位錯有限斷層模型[15],如圖4所示。

圖4 瑪多地震主震斷層同震位錯分布[15]Fig.4 Coseismic dislocation distribution of Maduo earthquake[15]

本文建立了2個有限元模型：均勻模型和分層模型,其材料參數如表2所示,楊氏模量和泊松比可由表中的密度、P波速度和S波速度計算得到。由于有限元模型的尺寸大大超過主震斷層,本文假定除模型上表面(地表)外,其余邊界上的法向位移為0。

表2 瑪多地震同震變形三維有限元分層模型的材料參數(自USGS)Table 2 Material parameters of the 3D finite element layered model of Maduo earthquake(USGS)

2.2 瑪多地震地表同震變形和同震應力變化

根據上述三維彈性有限元模型,本文模擬得到的瑪多地震地表同震位移場見圖5。東西向位移u大于零的部分位于斷層南側,小于零的部分主要集中于斷層北側,說明斷層南側以向東運動為主,斷層北側則以向西運動為主(圖5(a));南北向位移v大于零的部分位于斷層北側,小于零的部分主要集中于斷層東南部,說明斷層南側以向南運動為主,斷層北側以向北運動為主,呈現拉張的運動方式(圖5(b))。以上同震位移場分布特征與瑪多地震具有明顯的左旋走滑機制一致;垂向位移w大于零的部分位于斷層北部,小于零的部分位于斷層南部,斷層北部為上升盤,南部為下降盤,上升盤為斷層下盤,下降盤為上盤(圖5(c)),表明發震斷層具有正斷分量。

圖5 瑪多地震的地表同震位移分量等值線圖Fig.5 Contour maps of surface coseismic displacement components of Maduo earthquake

本文計算結果中瑪多地震引起的研究區域內主要斷層10 km深度的應力變化,如圖6所示。圖中正應力變化大于零的區域表示拉應力增加,主要位于斷層南西側和北東側,與拉應力減小的區域交錯分布(圖6(a));剪應力變化在發震斷層兩側是小于零的,在斷層尖端向外是大于零的(圖6(b));庫侖應力變化在斷層兩側基本上是小于零的,庫侖應力變化大于零的區域以斷層兩端為中心(圖6(c)和6(d)中的白色曲線分別是庫侖應力變化的0和0.001 MPa的等值線),在斷層兩側呈蝴蝶狀對稱分布(圖6(c)、6(d))。9～11 km深度范圍內除主震斷層帶兩側大于1.0級的余震基本位于10 km深度水平面上平行于主震的庫侖應力變化大于零的區域(圖6(d))。

黑色空心圓圈表示9～11 km深度的M≥1.0的地震。圖6 瑪多地震引起的深度為10 km水平面與主震震源機制一致的微元面上的同震應力變化Fig.6 The coseismic stress changes on the plane consistent with the focal mechanism of the main earthquake of Maduo earthquake at 10 km depth

本文同時計算了研究區域內與瑪多地震發震斷層相鄰的5條斷層上的同震庫侖應力變化,如圖7,計算使用的有效摩擦系數為0.4[24]。結果表明,瑪多地震引起鄂拉山斷裂南段、昆中斷裂東段、東昆侖斷裂帶的大部分(除中部靠西的小部分外)、瑪多—甘德斷裂北西段、南東段及主震震中附近以及達日斷裂帶北西段庫侖應力增加,地震危險性增加,這些地方是需要密切關注、加強地震監測的區域。

F1：達日斷裂帶;F2：瑪多—甘德斷裂帶;F3：東昆侖斷裂帶;F4：昆中斷裂帶;F5：鄂拉山斷裂帶,斷層數據引自詹艷等[4]。五角星為2021年5月22日青?，敹郙S7.4地震的震中。圖7 瑪多地震震中周圍主要斷裂上的庫侖應力變化Fig.7 Coulomb stress change on the major faults around the epicenter of Maduo earthquake

3 討論

3.1 瑪多地震地表同震變形與大地測量觀測數據及與前人計算結果的對比

為驗證本文采用分裂節點有限元方法模擬實際大地震同震變形的有效性,本文設計了瑪多地震的2種不同介質模型：均勻介質模型和分層介質模型,2種模型材料參數如表2。本文分別計算了2種介質模型的同震變形場,并將其處理為衛星視線方向LOS位移量,與GPS觀測的同震水平位移、InSAR升降軌LOS位移量觀測結果進行對比,結果如圖8所示。圖8(a)、8(b)中的黑線表示瑪多地震發震斷層,藍線表示有限元模型模擬結果,紅線表示GPS觀測結果(低頻),粉線表示GPS觀測結果(高頻),左側是均勻有限元模型模擬結果,右側是分層有限元模型模擬結果。圖8(c)、8(e)、8(g)分別是InSAR升軌的觀測數據、均勻和分層有限元模型模擬結果。圖8(d)、8(f)、8(h)分別是InSAR降軌的觀測數據、均勻和分層有限元模型模擬結果,表明均勻有限元和分層有限元模型得到的升軌和降軌的衛星視線方向LOS位移量與InSAR觀測得到的衛星視線方向LOS位移量有較好的一致性。

本文采用統計均方根(root mean square, RMS)方法對有限元模擬結果和觀測變形量之間的差異進行定量說明。RMS[15]定義如下

(8)

表3 瑪多地震同震位移的有限元模擬結果與大地測量觀測位移量之間的RMSTable 3 RMS comparison of finite element simulation results and geodetic observed coseismic displacement of Maduo earthquake

從表3可以看出,有限元模擬同震位移的結果與大地測量觀測結果之間的RMS在合理范圍之內,與楊君妍等[15]的計算結果一致,說明本文所采用的有限元方法對大地震同震變形的模擬是有效的。

本文通過分裂節點法將位錯數據加入有限元模型計算的同震位移結果表明,發震斷層南側以向東、向南及向下運動為主,斷層北側則以向西、向北及向上運動為主,總體表現為有正斷分量的左旋走滑斷層。該計算結果與楊君妍等[15]使用球形分層地球模型位錯理論、談洪波等[19]使用均勻半空間位錯理論進行同震位錯計算的解析法或半解析法所得結果基本一致。

3.2 瑪多地震對巴顏喀拉地塊及周邊地震危險性的影響

本文瑪多地震同震庫侖應力計算結果表明,發震斷層的兩側沿走向方向上庫侖應力下降,而在斷層兩端及以兩端為中心呈蝴蝶狀發散的區域內,庫侖應力增加。這個現象表明瑪多地震對發震斷層及斷層周邊地殼產生了顯著的應力調整：瑪多地震部分釋放了發震斷層所積累的應力(發震斷層兩側的庫侖應力變化大部分都是顯著降低),也造成斷層兩端及周邊部分區域的地殼應力積累(圖6(c))。該結果與楊光遠等[16]利用Coulomb 3.3程序在同等條件(10 km深度處和主震震源機制一致的微元面)下的計算結果對比,同震庫侖應力分布基本一致。我們的模擬結果還表明,9～11 km深度范圍內除主震斷層帶兩側大于1.0級的余震基本位于10 km深度水平面上平行于主震的庫侖應力變化大于零的區域(圖6(d))。研究區域內鄂拉山斷裂、昆中斷裂、東昆侖斷裂、瑪多—甘德斷裂以及達日斷裂這5條主要斷裂的同震庫侖應力結果表明：瑪多地震引起鄂拉山斷裂南段,昆中斷裂東段,東昆侖斷裂帶的大部分、瑪多—甘德斷裂西北段和東南段以及達日斷裂帶北西段庫侖應力增加(圖7)。鄂拉山斷裂南段與昆中斷裂東段位于巴顏喀拉塊體北部的東昆侖—柴達木塊體,歷史上該地區曾發生過M5～M6及M6～M7地震,此次瑪多地震引起該地區庫侖應力增加,具有較高的地震危險性,應對其密切關注。東昆侖斷裂帶歷史上曾發生過如1937年MS7.5地震(阿拉克湖—托索湖段和下大武段)、1963年MS7.0地震(秀溝—阿拉克湖段)以及2001年MS7.1地震(庫賽湖段),是一條活躍斷裂帶[25]?，敹嗟卣鹗宫斍摺斍螙|南段和托索湖—下大武段東南段庫侖應力增加,這與Li等[20]的結果一致。這2個斷裂段在瑪多地震時產生應力積累,加之瑪沁—瑪曲段近年來未發生大地震,導致該地區地震危險性大大增加,應當加強地震監測、做好防震防災工作?，敹唷实聰嗔阎卸螏靵鰬ψ兓植驾^為復雜,庫侖應力顯著增加與顯著降低的區段交錯分布,可能與其復雜的斷層結構有關[4]?，敹唷实聰嗔巡糠謪^段庫侖應力顯著降低可能與瑪多地震釋放了其部分應力有關,但其兩端的庫侖應力顯著增加,地震危險性有所提高。位于瑪多地震發震斷層南部的達日斷裂帶北段及中段庫侖應力增加,而其余部分庫侖應力水平下降,考慮到1947年達日M7.7地震所帶來的達日斷裂帶西北段庫侖應力增加及其較高的應力積累速率[26],其地震危險性仍需進一步研究與評估。

4 結論

本文采用分裂節點技術開發三維并行彈性有限元程序,并通過地震位錯模型的理論解驗證了有限元程序的正確性。該有限元程序能夠考慮復雜幾何、非均勻材料和位錯非均勻分布,可以模擬計算大地震同震變形。利用該程序模擬計算了瑪多MS7.4地震的同震位移場、應力場及庫侖應力變化,計算得到的地表同震位移結果與GPS和InSAR大地測量數據一致。模擬結果確認了瑪多MS7.4地震的發震斷層為一條左旋走滑兼有正斷分量的斷層。計算得到的庫侖應力增加的區域與瑪多地震的余震分布有較好的對應性。庫侖應力結果表明,瑪多地震引起了鄂拉山斷裂南段,昆中斷裂東段,東昆侖斷裂帶的大部分區段、瑪多—甘德斷裂西北段和東南段以及達日斷裂帶北西段庫侖應力增加,這些區域的地震危險性增加,是未來需重點關注的地區。

石耀霖院士和周永發博士對本研究自主研發的三維有限元程序編制給予了很多建設性的指導和幫助,周元澤教授對論文初稿提出了很多建設性的意見,在此一并表示感謝。