2022年門源 MS6.9 地震前斷層活動及應力狀態的數值模擬

李 媛 楊周勝 龐亞瑾 梁洪寶 劉 峽

1)中國地震局地質研究所,地震動力學國家重點實驗室,北京 100029 2)中國地震局第一監測中心,天津 300180 3)云南省地震局,昆明 650000

0 引言

據中國地震臺網測定,2022年1月8日1時45分,青海省海北藏族自治州門源縣發生MS6.9 地震,震中位于(37.77°N,101.26°E),震源深度10km,地表破裂約22km。此次地震的震中距1986年和2016年2次門源MS6.4 地震震中均約33km。3次門源MS>6.0強震的發生時間、空間距離均較近,發震斷層均為冷龍嶺斷裂西段。與前2次逆沖型地震(姜文亮等,2017)不同,本次地震幾乎為純左旋走滑型,與現有研究顯示的發震斷層的運動性質一致。開展發震斷層活動和應力狀態研究有助于促進對地震孕震機理的認識。

針對多次門源地震,一些學者在同震、震后形變特征及應力觸發等方面進行了研究(劉洋等,2019; 李振洪等,2022; 朱琳等,2022)。庫侖應力結果顯示,1986年門源MS6.4 地震及1927年古浪MS8.0 地震對2016年門源MS6.4 地震產生了應力加載效應(雷東寧等,2018); 2016年門源地震對2022年門源地震的發生有一定促進作用(李振洪等,2022; 朱琳等,2022)。但多期GPS應變結果(朱爽等,2022)顯示,2016年門源MS6.4 地震的影響主要集中在震中斷層附近非常小的區域,對2022年門源MS6.9 地震影響較小。本文認為,該局部斷層及附近地震如此頻繁,應力釋放遠高于其他區域,除與地震觸發有關外,應該還與發震斷層本身的運動狀態、長期應力累積趨勢和區域動力學背景有著必然聯系。

1986—2022年門源系列地震的發震斷層——冷龍嶺斷裂位于祁連山構造區。該區地處急劇隆升的青藏高原活動塊體向內陸擴展的前緣部位,受到青藏高原塊體長期的NE向推擠作用(鄧起東等,2002),加之受到剛性的塔里木和阿拉善地塊阻擋,新構造活動和變形強烈,斷裂構造十分發育。區內大型斷裂主要包括托萊山斷裂、冷龍嶺斷裂、昌馬-俄博斷裂、祁連山北緣斷裂、武威-天祝-莊浪河斷裂等,主要展布方向為NWW向和NW向,主要運動方式為走滑、逆沖或兩者的過渡類型(雷東寧等,2018)。研究表明,該區域的新構造活動以水平運動為主、垂直隆升運動為輔,兩者相差1個數量級; 水平運動和變形呈明顯NE向擠壓縮短、順時針旋轉和向SEE擠出的構造變形特征(Wangetal.,2020; 郝明,2012,王雙緒等,2013)。

關于冷龍嶺斷裂帶晚更新世以來的活動特征,早年的地質研究曾估計其滑動速率約為19mm/a(Gaudemeretal.,1995; Lasserreetal.,2002),而近年來的研究已將其量值限定在3.35～6.4mm/a的范圍內(何文貴等,2000,2010; 郭鵬等,2017)。在大地測量方面,基于高精度、長時間尺度的GPS觀測數據,多位學者(甘衛軍等,2005; Meade,2007; 李煜航等,2015)利用負位錯模型、三維彈性球面塊體模型和塊體-位錯模型反演了祁連塊體邊界主要活動斷裂的長期滑動速率、虧損滑動分布及閉鎖耦合變化,發現金強河斷裂、毛毛山斷裂和冷龍嶺斷裂的閉鎖程度高。一些學者基于三維有限元黏彈性模型,利用單期GPS數據,在青藏高原東北緣及海原-六盤山地區開展了斷層應力場狀態及與地震活動性相關的數值模擬研究(楊興悅等,2016; 石富強等,2018; 孫云強等,2018; 蔣鋒云等,2021),結果顯示,各個斷層的應力并不均勻,強震主要發生在高應力狀態的斷層上,當區域應力累積到較高水平時,可在短期內經歷一系列地震或地震叢集?，F今祁連山構造區的木里-江倉斷裂、金強河-老虎山-毛毛山斷裂的剪應力積累較大。龐亞瑾等(2017)利用多期GPS數據,基于三維彈性模型研究了九寨溝MS7.0 地震前不同時段的構造應力積累情況,結果顯示,現今的地殼應力積累過程有利于區域內左旋走滑型地震的發生。前人的研究為了解該區域地震過程中斷層的運動、應力狀態及變化特征提供了寶貴的資料和經驗。但是,在確定斷層不同區段的滑動速率方面,地質研究數據仍有待補充。由于在以往的數值模型中簡化了斷層形態,且受到研究側重點不同等因素的制約,目前針對冷龍嶺斷裂在不同時間尺度上的斷層運動變化及其與地震關系的研究還較少。鑒于模型的精細度將很大程度影響計算效率和結果,考慮到門源地震的發震斷層在走向展布上的特殊幾何形態,如果能提高模型的精細程度,對提升地震孕震過程中斷層運動和應力特征的研究認識將會有較大的幫助。

基于前人工作的基礎,本文利用地震前1991—2015期和2017—2021期的GPS速度場,綜合考慮各級地質單元、地殼的分層結構、介質的橫縱向不均勻性、復雜的斷層系統及斷層形態等多種因素,建立了祁連山構造區精細的三維有限元黏彈性模型,計算分析了祁連山構造區在長期的構造運動環境下區域應力積累的基本格局,區域內斷層的長期滑動速率、應力累積速率,以及這些量值在門源MS6.9 地震前約5a時段內的變化特征。同時,結合跨斷層短水準、震源機制解等資料,探討門源MS6.9 地震前近震區強震發生的動力學背景和孕震環境。

圖1 祁連山斷裂帶構造背景圖

1 有限元數值模型構建

1.1 區域及活動斷裂帶的三維模型

為充分利用研究區內的GPS觀測數據約束模型,同時考慮塊體運動差異對斷層運動的影響和模型計算的邊界效應等,將模型擴大為青藏高原東北部的整個區域。如圖2a 所示,三維模型考慮了地殼的分層結構、介質的橫縱向不均勻性因素,涵蓋了多個活動地塊,東至鄂爾多斯地塊腹地,東南至四川盆地,西北至阿拉善地塊腹地,中南部為祁連山地塊。各地塊邊界以中國大陸活動地塊劃分方案(張培震等,2003)為依據。模型厚度為100km,分為上、中、下地殼和地幔層,地殼分層結構主要依據全球地殼模型Crust1.0(Laskeetal.,2013)。由于莫霍面起伏對模型結果具有顯著影響,且Crust1.0數據的精度較低(1°×1°),模型的殼幔分層界面采用精度較高的中國海陸1︰500萬莫霍面深度數據(郝天珧等,2014),地殼厚度自西南的祁連山地塊向東北的阿拉善-鄂爾多斯地塊逐步遞減。

圖2 三維有限元模型的邊界條件和擬合殘差結果

同時,模型考慮了復雜的斷層系統,包含30余條大型斷裂,其中祁連山斷裂帶主要包括黨河南山斷裂、陰凹槽斷裂、托萊山斷裂、冷龍嶺斷裂、昌馬-俄博斷裂、祁連山北緣斷裂、榆木山斷裂、龍首山斷裂、金強河斷裂、武威-天祝斷裂、莊浪河斷裂等。斷層形態主要依據中國斷層圖(鄧起東,2007)詳細刻畫。如表1所示,大部分斷層的傾角為60°～80°或更為陡立,以走滑運動為主。參考前人的研究(楊興悅等,2016; 龐亞瑾等,2017; 石富強等,2018),不考慮斷層介質橫向不均勻性的影響,同時為兼顧計算收斂問題,將斷層統一簡化為寬約3km的直立弱化帶,埋深至上地殼底部深度。

表1 模型物性參數設置

模型采用ANSYS軟件中具有 2 次位移模式的Solid187六面體單元進行網格劃分。斷層單元的邊長尺寸為3km,其他單元控制在5km以內。模型共有379633個單元,71080個節點。網格劃分后的三維模型見圖2a,其中圖框部分為本文重點討論的祁連山構造區。規定模型的xyz坐標系為:x軸向E,y軸向N,z軸向上。

1.2 模型物性參數的選取

模型中,上地殼采用脆性結構,其他分層采用黏彈性體,介質參數包括楊氏模量、泊松比和黏滯系數。在假定介質彈性各向同性的前提下,基于波速結構結果(Stolketal.,2013)和Crust1.0(Laskeetal.,2013)巖石密度結果,依據式(1)(張東寧等,1989)計算不同地塊的楊氏模量和泊松比。結果顯示,各塊體間的彈性性質在橫向上的差距較小,與前人的計算結果相符(Xiaoetal.,2015; 楊興悅等,2016; 祝愛玉等,2019)。各地塊在上地殼、中地殼、下地殼及上地幔的楊氏模量、泊松比和密度設置如表1所示。參考前人對斷層的設置方法,將斷層的楊氏模量設置為較周圍介質偏小約1個量級,即1×1010Pa,泊松比較周圍介質高0.03(Xiaoetal.,2015; 楊興悅等,2016; 祝愛玉等,2019)。關于介質黏滯系數,相關研究表明,模型中不同地塊的巖石圈流變性質在橫向和縱向上均有差異,且在中國西部低黏滯系數的下地殼普遍存在,下地殼的黏滯系數比巖石圈地幔一般低1～2個數量級,尤其是青藏高原存在著強度最低的下地殼(石耀霖等,2008; 陳連旺等,2011; 孫玉軍等,2013)。本文采用祝愛玉等(2019)和孫玉軍等(2013)給出的黏滯系數結果。

(1)

其中,ρ為巖石密度,VP、VS為橫波速度及縱波速度,?為泊松比,E為楊氏模量。

1.3 模型加載和計算

為探索不同時間尺度形變特征,本文采用門源MS6.9 地震前1991—2015期(Zhengetal.,2017)和2017—2021期的GPS速度場作為模型邊界加載條件。其中,2017—2021期GPS水平運動速率場的數據處理和解算由最新的GAMIT/GLOBK軟件(2)http:∥geoweb.mit.edu/gg/。完成。在該時間段內青藏高原東北緣發生了九寨溝MS7.0 地震,在速度場中去掉了該地震的同震影響。

選取位于模型邊界附近的GPS測站,將其觀測結果通過雙次樣條插值至模型側面作為水平向運動約束。注意在GPS站點的選取過程中,需要剔除與周邊測點運動顯著不協調的測點,以降低模型邊界附近的斷層活動等局部運動信息對模擬運算的影響。參考前人研究(劉峽等,2010; 石富強等,2018; 萬永魁等,2021),設定水平速度不隨深度改變。另外,模型的上表面為自由面,下表面水平向自由、垂向固定。

本文采用靜力學分析,巖石圈視為線彈性體和理想黏彈性體串聯的三維黏彈性模型,其變形服從Maxwell體的本構關系(Liuetal.,2003)。Maxwell體在短時間內表現為彈性體,隨著外力作用時間的增加,儲存在彈性體內的勢能逐漸消失在黏性體中,隨后應力松弛。弛豫時間一般定義為黏滯系數與楊氏模量的比值。當流變時間遠大于弛豫時間時,應力得到松弛并趨于穩定(Caoetal.,2009)。本文模型的弛豫時間范圍為114～1655a。本文將加載時間設定為30000a。第1步,施加1991—2015期GPS速度場約束,開始采用的步長為100a,持續加載至29800a后將步長改為1a,繼續加載195a,此步計算的目的是得到一個反映較長時間尺度、相對穩定的運動場和應力場; 第2步,在第1步的基礎上繼續施加2017—2021期GPS速度場約束,加載時間為5a,步長為1a,其計算目的是獲取2022年門源MS6.9 地震發生之前約5a內研究區運動場和應力場的變化。在計算過程中可以觀察到,當加載至約10000a時,模型地殼和地幔應力已趨于平穩。計算采用Newton Raphson迭代算法進行求解,將計算劃分為一系列載荷增量步,在每個步長內自動劃分出多個子步長; 采用殘余力和位移的范數控制收斂,收斂精度為0.1%。需要說明的是,本文計算顯示,長期穩態構造加載下的應力場平均累積速率約為1kPa/a。而該區域歷史大地震對區域應力加卸載的影響量級在10～100kPa(萬永革等,2007; 鄧園浩等,2018; 湯大委等,2023),相比于萬年尺度的構造應力場而言,量值非常小,故本文在加載中忽略地震事件的影響。

2 數值模擬結果分析

2.1 數值模擬計算的合理性分析

圖2b為1991—2015、2017—20212期GPS觀測值和模型表面邊界的加載。如圖2c 所示,模型范圍內198個GPS測點的觀測值與模擬值的殘差方向較為散亂,呈無序狀態,且殘差值大小不一,說明其中不包含顯著的整體運動信息。對模型內GPS測點的N向和E向速率殘差值進行統計分析,結果顯示(圖2d,e),2期殘差值基本服從正態分布,殘差<2mm的測站分別約占總測站數的90%和81%,其中殘差<1mm的測站分別占比約76%和68%。綜上,本文認為2期的邊界約束和模擬結果比較合理。

此外,表2 給出了本文第1載荷步計算獲取的各斷層的運動性質與速率的模擬結果,作為對比,也列出了前人的地質考察結果。模擬得到的斷層性質與前人的地質研究結果相符,模擬給出的大部分斷層的滑動速率與地質結果相近,進一步印證了模擬計算的合理性。

表2 祁連山構造區主要活動斷裂幾何形態和活動特征表

2.2 1991—2015期區域應力場模擬結果揭示強震在長時間尺度上的孕育背景

應力場可直觀地反映地殼受力狀態及地震孕育的背景。圖3a 為利用較長時間跨度的1991—2015期GPS速度約束,計算給出的研究區上地殼10～15km范圍的最大、最小主應力和最大剪應力累積率結果,為應力的平均年變化率,代表了構造應力場長期的變化趨勢。根據該圖,應力場自西向東逐漸由NNE向擠壓和NWW向拉張,順時針旋轉為NE向擠壓和NW向拉張,主壓應力方向大致與斷層走向垂直,與中國地殼應力場(謝富仁等,2004)中的最大主應力方向相符。此外,本文收集整理了祁連山地區1970年1月—2021年9月利用HASH和CAP方法進行解算和精定位的,共619個3級以上地震震源機制解結果,基于斷層面解的阻尼應力張量反演法(Hardebecketal.,2006),獲得水平構造應力場,見圖3b。本文的應力模擬結果為上地殼網格單元的應力結果,深度范圍為10～15km,而震源機制解涉及的深度則更寬泛,盡管2種應力結果在某些局部存在差別,但總體上大致符合。值得關注的是,模擬結果顯示在門源MS6.9 震源區及附近,應力場處于NE向擠壓和NW向剪切共同作用的動力學背景之下,為主壓應力場由NNE向向NWW向轉換的過渡地區,與GPS應變率場研究結果(Zhengetal.,2013; Wangetal.,2020)一致?？梢?現今應力場的長期的、整體的積累趨勢,促進了NWW向斷層的走滑和擠壓運動,有助于托萊山斷裂、冷龍嶺斷裂等左旋走滑斷層上走滑型和逆沖型地震的孕育。

圖3 研究區上地殼主應力累積率場分布

圖3a顯示,研究區自NW向SE最大剪切應力累積率逐漸減小,西部地區最大達2.4kPa/a,東南地區為0.6～1.2kPa/a。GPS研究給出的1999—2017年最大剪應變率結果(朱爽等,2022)同樣顯示了西強東弱的特征,與本文剪應力累積率相符合。此外,在托萊山和冷龍嶺斷裂交會部位,即門源MS6.9 震源區附近的剪應力累積率較周圍大,最大為2.0kPa/a,其東西兩側約為1kPa/a,揭示門源MS6.9 震源區受到強烈的剪切作用,更容易發生地震。

2.3 1991—2015期斷層長期平均運動與應力場模擬結果

本文模型利用連續、具有一定寬度、黏滯系數和楊氏模量都顯著低于周圍穩定地塊的薄層單元,模擬斷層的運動學特征和釋放能量的力學性質,與實際斷層的非連續錯動滑移仍有區別。盡管如此,從模擬結果看,斷層兩側的位移仍然有顯著差異(圖4)。參考已有的處理方法(祝愛玉等,2015),本文的斷層擠壓和走滑速率分量分別是斷層兩側網格節點垂直于斷層走向和沿斷層走向的相對位移速率,相應地,斷層應力累積率為垂直于斷層走向和沿斷層走向的法向正(擠壓)應力和剪切應力累積率。

圖4 斷層長期滑動分量與應力積累速率模擬結果

圖4a、b為利用較長時間跨度的1991—2015期GPS速度場作為約束,計算給出的斷層長期平均滑動速率。結果顯示,區內斷裂活動總體表現為:大部分NWW向斷裂受到左旋走滑兼擠壓作用,走滑速率為0～5mm/a,擠壓速率為0～3mm/a; NNW向斷裂表現為右旋走滑兼擠壓,如莊浪河斷裂和榆木山斷裂東段,與地質考察結果較為符合(國家地震局地質研究所等,1993; 侯康明等,1999; 袁道陽等,2002)。斷層應力累積率結果(圖4c,d)顯示,力的性質與運動性質相匹配,大部分NWW向斷層既受到擠壓正應力又受到左旋剪切應力,NNW向斷層受右旋剪應力。上述結果與楊興悅等(2016)、戴華光等(1996)給出的有限元結果及李煜航等(2015)根據GPS計算的主要大斷層的滑動速率一致,與該區域晚新生代以來處于長期擠壓與剪切共同作用的動力學背景相符。

空間上,斷層運動和應力呈現不均勻分布。研究區內斷層的運動趨勢為:自西向東擠壓速率逐漸減小,由3.0mm/a減小為0.2mm/a; 自西向東走滑速率逐漸增大,由0.9mm/a增加為4.6mm/a。這揭示了祁連山構造帶在青藏高原東北緣構造變形中的轉換和調整作用,且尤其以托萊山斷裂—冷龍嶺斷裂的斷層運動和應力分段性特征最為明顯:從托萊山斷裂西段至冷龍嶺斷裂東段,擠壓速率從2.5mm/a逐漸減少為0.2mm/a,而走滑速率從1mm/a逐漸增大至4.6mm/a; 剪切應力積累年變化率為0.5～3.8kPa/a,擠壓正應力積累年變化率為0.2～1.0kPa/a。從圖4b 中可以看出,相對周邊斷層,托萊山斷裂和冷龍嶺斷裂左旋走滑運動較強,但在冷龍嶺斷裂西段的幾何拐折部位,擠壓運動分量和走滑分量較兩側斷層段落更弱,而斷層的法向擠壓應力和左旋剪應力積累率較大,表明運動受阻。

2.4 2017—2021年模擬斷層運動與應力增量場的結果分析

在2017—2021時段GPS速度場的約束下,計算得到了2022年門源MS6.9 地震發生之前約5a內的斷層滑動速率和應力累積率分布。為了更加直觀地反映斷層運動的變化情況,將2017—2021期的模擬結果與1991—2015期結果相減,得到斷層滑動速率和應力累積率的變化量,即增量(圖5)。整體來看,擠壓、走滑速率增量的絕對值范圍分別為0～1.5mm/a和0～2.5mm/a。大部分斷層段的變化量并不大,在1mm范圍內,局部斷層段變化較大,如榆木山斷裂東段、黨河南山斷裂中段和冷龍嶺斷裂西段。斷層的增量分布表現為顯著的分段差異和不均勻性:一方面,走滑速率增量(圖5b)顯示,以101°E為界,西部大多數斷層的走滑速率增量為正值,東部大部分斷層為負值,表明在2022年門源MS6.9 地震前,西部斷層左旋走滑運動有所增強,而東部則有所減弱,尤其是冷龍嶺斷裂西段減弱較為顯著,量值約為1mm/a,而此段落的剪應力積累速率增量(圖5d)為正值,即斷層的剪切應力積累加速; 另一方面,圖5a,c顯示,擠壓速率和正應力累積率(擠壓速率和擠壓應力為負值)增量為正的斷層分段主要包括黨河南山斷裂、陰凹槽斷裂、托萊山斷裂東段、冷龍嶺斷裂、祁連山北緣斷裂中段及NNW走向的斷層,表明近期這些斷層的擠壓運動有所減弱,壓應力積累減緩。尤其在托萊山斷層東段及冷龍嶺斷層西段鄰接的部位,近期斷層的滑動速率和壓應力累積率顯著減弱,這一變化有利于斷層錯動。

圖5 震前約5a較長期的斷層滑動速率與應力累積率變化結果

3 討論

3.1 發震斷層長趨勢運動和應力累積特點與門源 MS6.9 地震的關系

根據本文利用1991—2015期GPS結果給出的區域應力場、斷層長期的運動和應力累積分布,整體而言,托萊山和冷龍嶺斷裂受剪切作用力較強(圖3),斷層左旋走滑速率達4.6mm/a,與野外地質考察顯示的這2條走滑型斷裂具有4.0～5.38mm/a和3.36～4.62mm/a的水平滑動速率(何文貴等,2000,2010)差別不大。但在托萊山斷裂和冷龍嶺斷裂鄰接幾何拐折部位及附近區域,擠壓速率和走滑速率比相鄰的斷層分段更弱(圖4a,b),同時該部位的斷層法向擠壓應力累積率和左旋剪應力累積率較大(圖4c,d),這與該斷層段可能處于閉鎖狀態(葉茂盛等,2018; Guoetal.,2022; 劉雷等,2023)有關。在閉鎖階段,斷層運動速率很小而能量快速積累,直到發生一次或若干次較快滑動,將積累的能量迅速釋放掉?？梢酝茢?這樣的斷層區段更容易發生地震。

中國大陸板內地震是在統一的、大的構造應力場作用下發生的(Tapponnieretal.,1977),主要受發震斷層段的構造幾何形態、介質物性、斷裂活動習性和應力積累程度等孕震環境控制(陶瑋等,2011; 張培震等,2013)。鑒于祁連山地塊內部地殼結構和介質物性參數橫向差異不大(Laskeetal.,2013; 郝天珧等,2014),推斷引起冷龍嶺斷裂西段滑動和應力積累減弱的因素主要與該斷層分段走向所呈現的特殊形態緊密相關。過去的研究多關注塊體邊界斷裂,且大多對其形態進行簡化。而1986—2022年門源系列地震發生在冷龍嶺斷裂西段走向呈現拐折形態的區域,具有特殊性,采用簡化的斷層模型難以解釋其孕震過程的特征性變化。本文構建的連續三維黏彈性有限元模型(圖2a 所示)考慮了各級地質單元、殼幔分層、地殼的橫縱向不均勻性和介質的黏彈性等多種因素,同時構建了復雜的斷層系統,細化了斷層的幾何形態,所獲得的區域形變場、應力場和斷層運動,應是周邊塊體相對運動、巖石層的黏彈結構和斷層幾何形態和空間展布等共同作用的結果。與已有模型的結果相比(楊興悅等,2016; 石富強等,2018; 孫云強等,2018),本文模擬的斷層運動及其應力分布具有顯著的分段性特征,體現了斷層的幾何形態對斷層活動的重要影響。

除2022年門源MS6.9 地震外,許多強震發生在斷層拐折附近(Andrews,1994)。例如,San Andreas大斷裂分別于1690年、1630年和1857年發生了3次MS>7.0地震,震中均位于斷層拐折附近; 1976年唐山MS7.8、2013年蘆山MS7.0 及1927年古浪MS8.0 地震同樣為發生在拐折斷層上的事件(Schwartzetal.,1990; Aochietal.,2002; 張竹琪等,2013)。Andrews(1994)指出斷層的拐折部位對斷層的滑動具有抑制作用; 薛霆虓等(2009)研究表明,幾何彎曲的斷層導致了應力的集中,斷層拐折附近剪應力顯著增加且多發育強震。針對拐折斷層強震多發的位置,通過物理實驗和數值模擬研究表明:地震并未發生在拐折顯著的位置,而是發生在拐折附近或拐折兩側的平直斷層處,且拐點起到類似閥門的作用(馬瑾等,2012)。而1986—2022年3次門源MS>6.0地震的震中正位于拐折斷層的兩端靠近平直斷層分段的位置。冷龍嶺斷裂帶西段特殊的形態,很大程度上導致了該局部斷層長趨勢的運動虧損,與兩側的斷層段運動存在顯著不協調,促使應力在該斷層段快速積累。因此,冷龍嶺斷裂帶西段具有強震多發的孕震背景和構造條件。

3.2 發震斷層震前中短期運動和受力特點與2022年門源 MS6.9 地震的關系

圖6 為2022年門源MS6.9 地震前,祁連山構造區跨斷層短水準觀測的時序曲線,所用數據經過去周期和4～6階多項式趨勢擬合等處理。曲線趨勢向下代表所跨斷層出現擠壓、逆沖運動,向上則為拉張、正斷運動。將跨斷層短水準表現出的張壓運動特征,與2017—2021期相對于1991—2015期擠壓速率增量的模擬結果(圖5a)進行對比,發現兩者具有較為一致的變化趨勢。例如,祁連山北緣斷裂東段上的扁都口、石灰窯口和冷龍嶺斷裂東段的白土莊跨斷層場地短水準測段在2017年后均表現為擠壓運動偏離線性和趨勢減緩的特征; 大野口測段2018年以來為拉張增強或擠壓減弱的運動趨勢,與此同時,模擬結果(圖5a)顯示這些跨斷層場地所處的斷層分段,其擠壓速率增量均為正值,表現為與跨斷層水準曲線相符的擠壓減弱變化。再如,金強河斷裂上的烏鞘嶺、松山跨斷層短水準在2018年以來呈擠壓增強的運動趨勢,模擬結果(圖5a)也表明該斷層的擠壓速率增量為負,近期擠壓運動顯著增加。又如,昌馬-俄博斷裂東段的童子壩和龍首山斷裂東段的滑石口井短水準自2009年以來擠壓運動速率幾乎未發生改變,而模擬結果(圖5a)顯示相應斷層分段的擠壓速率增量接近零值。上述跨斷層形變實測數據可與模型計算結果相互印證。

圖6 震中距150km以內的異常測段時序曲線

2017—2021期相對于1991—2015期的斷層運動速率增量結果(圖5a,b)顯示,2022年門源MS6.9 震前中短期,近震區斷層的擠壓運動和走滑運動有所減緩,尤其是托萊山斷裂東段和冷龍嶺斷裂帶西段,其走滑速率由原來的2.0mm/a減少至1.0mm/a,擠壓速率由原來較弱的運動速率1.5mm/a減少至0.3mm/a,說明斷層擠壓和剪切變形在震前中短期進一步受阻,可能與地震前冷龍嶺斷裂存在閉鎖現象(Guoetal.,2022)有關。利用GPS剖面分析得到的斷層滑動速率和區域應變場動態結果(朱爽等,2022)也顯示,冷龍嶺斷裂帶的左旋走滑速率和區域剪切應變率在2017年后有所減弱。震前中短期運動趨勢有減緩的現象,在1976年唐山MS7.8、1989年美國LomaprietaMS7.1、2013年岷縣漳縣MS6.6、2008年汶川MS8.0 等地震發生前也曾出現(Lindeetal.,1992; 杜方等,2009; 牛安福等,2020)。此外,2017—2021期相對于1991—2015期斷層剪應力累積率的增量結果(圖5d)顯示,冷龍嶺斷裂西段的剪應力積累速率進一步增強,擠壓正應力積累速率有所減緩,剪應力、正應力的變化有利于左旋走滑型地震的發生。上述現象說明震源區在震前中長期已處于高應力積累水平,隨著地震臨近,發震斷層面法向擠壓應力緩慢增高,剪切摩擦力同步增大,導致斷層運動“遲滯”(李媛等,2023),促使運動趨勢偏離線性和減緩。

不僅如此,2017—2021期相對于1991—2015期的增量結果也顯示,在與冷龍嶺斷裂西段彎曲構造相似的六盤山東麓斷裂的彎曲部位,斷層擠壓和走滑運動剪切變形在近期略有減少,運動進一步受阻,可能預示著六盤山東麓斷裂的彎曲部位附近的強震危險性有所增加。

4 結論

本文利用1991—2015期和2017—2021期GPS速度場,通過建立青藏高原東北部三維有限元模型,計算了地殼形變和應力積累的長趨勢特征及其在門源震前中短期的變化。結合震源機制解、跨斷層短水準等結果,綜合分析了門源MS6.9 地震近震區強震多發的孕震背景和震前形變特征,對其機理進行討論。得到如下結論:

(1)在青藏高原NE向的推擠作用下,現今門源地震近震區受到NE-SW向擠壓和NW-SE向拉張的構造應力場作用,有助于促進該區域內NWW向斷層的左旋走滑和擠壓運動,利于左旋走滑斷層上地震的孕育。

(2)冷龍嶺斷裂西段具有強震多發的孕震背景和構造條件:1)門源MS6.9 地震震中局部區域最大剪切應力累積率比周圍區域更大; 2)冷龍嶺斷裂西段呈幾何拐折形態,其特殊的構造導致不同的斷層分段之間應力積累的差異性和運動不協調,與兩側斷層分段相比,該幾何拐折部位的應力積累速率較高,但走滑運動在空間上受到了一定阻礙,斷層速率偏低,運動存在虧損; 3)自2017年以來至2022年門源MS6.9 震前,冷龍嶺西段斷層剪切應力積累速率加快,斷層運動尤其是走滑運動顯著減弱,斷層運動進一步“遲滯”。上述斷層運動和應力狀態有利于冷龍嶺斷裂西段上地震的孕育和發生。綜上所述,冷龍嶺斷裂西段具有高應力積累狀態和斷層運動虧損的孕震背景,判斷該斷裂分段未來依舊是強震頻發的危險區段。

此外需要指出的是,受計算量、收斂問題等多方面的限制,本文模型僅是在幾何形態(即走向)上對斷層進行了詳細處理,不考慮傾角,將斷層簡化為直立斷層; 同時,模型邊界約束沿殼幔深度方向并無變化。這可能導致斷層面傾滑運動分量的模擬值偏小,因此仍有待進一步改進。

致謝跨斷層數據來源于中國地震局第二監測中心和甘肅省地震局; 本研究使用的有限元軟件為中國地震局地質研究所提供的正版Ansys商用有限元軟件; 馬勝利研究員對文章給予了多次指導; 審稿專家為本文提出了寶貴建議,使得文章更加完善。在此一并表示感謝!