厄瓜多爾俯沖區震間閉鎖與粘彈性變形研究

顏復康 田 鎮 楊志強 楊 兵 梁 沛

1 長安大學地質工程與測繪學院,西安市雁塔路126號,710054

板塊俯沖是地殼向地幔進行物質和能量傳輸的重要方式之一,也是地殼縮短、消亡的主要因素[1-2]。長期的板塊俯沖作用會沿著俯沖界面產生應力積累,并在震空區逐漸增加,當應力超過巖層承受強度時,則會以地震、海嘯、火山噴發等形式釋放,給人類生命財產安全帶來損失[3]。

厘清巨型逆沖斷裂閉鎖特征與地表形變的關系,是理解俯沖帶地震周期模式的關鍵環節[4]。根據現代大地測量資料,反演斷裂面閉鎖程度和滑動虧損分布特征是評估斷裂帶未來地震危險性的重要手段[5]。研究表明,在震間深部蠕滑區,俯沖板塊的推擠作用會引發軟流層緩慢的粘性流動,從而導致地表速度場發生顯著變化,并展現出更遠距離的長波形變模式,即震間深部粘彈性圈層的動力效應,會在一定程度上影響地表變形[4,6-8]。然而,目前基于Okada彈性半空間的位錯理論在推估斷層閉鎖特征時往往會忽略板塊粘彈性變形隨時間的變化特征,存在一定局限性[9-11]。同時為擬合GPS實測形變場,彈性模型可能會在一定程度上將斷層的閉鎖區域估算得更深,從而導致區域閉鎖特征估計錯誤[8,12]。因此,基于考慮深部介質差異的粘彈性模型來重新評估俯沖區的閉鎖深度十分重要。

研究斷裂的震間閉鎖深度能夠有效判斷斷層應力積累特征,對于評估區域未來地震危險性和研究巖石圈流變特性都具有重要意義[13-16]。前人針對厄瓜多爾俯沖區的研究大多基于彈性回跳理論來反演斷層閉鎖特征,鮮有從粘彈性模型角度出發,并基于數值模擬的方法研究該區域的斷層閉鎖分布[17]。基于此,本文以該區較為密集的GPS速度場為約束,分別建立厄瓜多爾俯沖區純彈性與彈性+粘彈性的地球模型,并采用有限元數值模擬方法,研究深部粘彈性圈層對震間地表變形的響應特征,從而對俯沖區的閉鎖深度重新進行評估。

1 地質背景

南美板塊位于納斯卡板塊、可可斯板塊與大西洋板塊之間,受到來自納斯卡板塊東向強烈沖擊和大西洋海底拉張應力,導致板內形變非常復雜[13,18]。

本文研究區位于納斯卡板塊和南美板塊俯沖區北部,主要包括厄瓜多爾與哥倫比亞南部,也稱厄瓜多爾俯沖區,西側的納斯卡板塊以約55 mm/a的速度沿N83°E方向俯沖至厄瓜多爾的南美板塊之下[16]。大地測量數據表明,在厄瓜多爾北部,Rio-Chingual-La-Sofia斷層存在約7.6±0.5 mm/a的右旋運動;而中部Pallatanga斷層右旋運動減少到5.3±0.4 mm/a;南部則幾乎無應變積累,GNSS速度降至約3 mm/a,并轉向東南方向[1,17,19]。由此可知,厄瓜多爾地區地殼運動的主要動力來源于納斯卡板塊和南美板塊之間的俯沖碰撞。此外,在厄瓜多爾南部,位于納斯卡板塊之上的卡內基海脊與南美洲大陸發生次生碰撞,導致局部構造變形極為復雜[20]。為簡化模型,本文在之后研究中并未考慮卡內基海脊對南美大陸的碰撞變形作用。

2 數據和方法

2.1 GPS速度場

Mocquet等[11]采用GAMIT/GLOBK軟件和地球動力學的標準策略,獲取了1994～2012年GPS時間序列,在分析所有GPS站點速度時考慮GPS時間序列中與時間相關的噪聲。本文采用的GPS速度場來自該研究中相對穩定的南美洲大陸水平速度場(圖1),區域速度場圖像表明,南美板塊整體呈東向運動,并沿內陸方向運動速率逐漸降低。

由于卡內基海脊海溝對厄瓜多爾南部存在碰撞作用,導致鄰區GPS數據出現明顯的離群現象。為簡化模型,本文并未顧及海脊的作用機制,因此剔除該區域部分GPS站點,重點考慮納斯卡俯沖板塊對南美大陸的俯沖作用,最終選取的GPS站點如圖1所示。

2.2 有限元模型

為厘清震間地表形變與深部粘彈性圈層的響應機制,本文分別建立純彈性模型與彈性地殼+粘彈性軟流圈模型。其中,有限元模型主要考慮海洋板塊和大陸板塊,在板塊上方被定義為延伸至各自估計彈性厚度的純彈性體,分別為大陸地殼和海洋地殼;在板塊下方則認為巖石圈和軟流圈地幔為各向同性的大陸地幔域和海洋地幔域,這些地幔域在彈性模型中具有純彈性性質,而在粘彈性模型中則具有粘彈性性質[8]。本文建立的三維粘彈性有限元模型可分為4個塊體,分別為海洋地殼、陸地地殼、海洋地幔和陸地地幔(圖2(a)),模型幾何構建與地殼厚度主要來自crust1.0(https:∥igppweb.ucsd.edu/)與部分前期研究成果[10,17,21]。由于莫霍面的復雜性,本文對所建模型進行一定程度的簡化處理,將海洋地殼與陸地地殼厚度均設置為30 km,該深度以下為粘彈性地幔區域[8]。

圖2 研究區有限元模型示意圖

在建立的三維有限元模型中,俯沖界面數據從Slab 2.0(https:∥www.usgs.gov/)獲取。由于獲取的俯沖界面深度約為400～500 km,為降低數值模擬中邊界效應對研究區的影響,同時考慮到地幔深處幾何界面對模擬的影響較小,俯沖界面延伸至800 km深度[22]。同樣,對于研究區周邊區域也采用擴張方法,使整個有限元模型明顯大于本文研究區。圖2(b)為構建的研究區幾何模型,模型長度為2 600 km,寬度為1 100 km,深度為800 km。在本文模型中,所有橫向邊界均固定為在法線方向具有零位移,但在切線方向不受約束;底面邊界在垂直方向固定,但在水平方向不受約束;頂部邊界在水平和垂直方向均不受約束。該設置會使斷層沿板塊界面滑動及地幔粘彈性松弛產生的應變在板塊內部累積[23]。

在對三維實體模型進行格網劃分時,為更好地計算斷層破裂附近的位移精度、數值計算精度與收斂性等[24-25],突出重點研究區厄瓜多爾及海溝的構造變形,本文首先將斷層面上海溝附近區域劃分為邊長約2 km的四面體單元,設置格網尺寸隨距斷層面的距離以1∶1.2比例逐漸增大,最終共劃分734 072個四面體單元,模型共包含143 157個節點,最小格網尺寸為2 km,最大格網尺寸為110 km(圖2(b))。該方法可避免四面體劃分出現高度扭曲,從而出現非收斂解和數值誤差,同時在深地幔中使用粗元素可節省計算時間和成本,且不影響分辨率精度[23-24,26]。

本文忽略震間模型的重力因素,將其視為絕對應力狀態的擾動,這些應力擾動通過沿斷層界面以運動學方式指定反向滑移率進行模擬[8]。同時,為使模擬結果能較客觀地反映該區域構造變形特征,建模時忽略地表曲率[27]。有限元節點上的累積位移可根據其在假定笛卡爾坐標系內的原始位置進行計算。即使在遠場中,地球曲率導致的水平差異也可以忽略不計,球面坐標系和笛卡爾坐標系測試模型之間的差異很小[23,27]。此外,已有研究表明,在模擬百年尺度的震間變形時,該區域歷史大地震的影響較小,因此本文在之后研究中并未考慮歷史地震的影響[8]。

2.3 巖石圈流變結構配置

本文采用的P波、S波速度及地殼、地幔密度、粘彈性等參數均來自文獻[12]和[25]。彈性模型和粘彈性模型之間唯一的區別是地幔域是否施加粘性參數。在彈性模型中,地殼與地幔變形視為彈性回彈,無粘度系數;而在粘彈性模型中,地幔變形視為彈性變形和與時間相關的粘性行為的組合,大陸地幔與海洋地幔的粘度分別為4×1019Pa·s和1×1020Pa·s[10]。粘彈性材料的彈性特性(如楊氏模量、泊松比和密度)與對應彈性模型中的彈性材料相同[10]。表1為本文模型的地塊屬性。

表1 地塊屬性

對于彈性模型,由于彈性變形是與時間無關的變形行為,因此本文模擬時間尺度為1 a的地表位移。而對于具有隨時間變化的變形行為的粘彈性模型,本文使用自適應時間步長方法模擬200 a的震間閉鎖[8]。該方法基于本構模型和變形率返回穩定時間步長,從而捕獲粘彈性材料的穩定響應[8]。對于粘彈性模型,地表速度場是從模擬的最后一個時間步長計算而來,因此該時間之后的粘彈性響應(約為Maxwell弛豫時間的10倍,已經穩定)速度可以代表地震后期速度。

本文采用Pylith有限元計算軟件模擬厄瓜多爾震間閉鎖效應在不同模型下的最佳閉鎖深度[28]。對于大地震震后幾百年的震間變形來說,雙粘度Burgers材料所代表的瞬態流變學對實驗的影響微乎其微,甚至對總的粘彈性貢獻有所低估。而南美洲大陸板塊下方海洋納斯卡板塊快速俯沖引發的大型(MW>7.5)大推力地震,其特征重現時間為100～250 a[29]。因此,本文選擇忽略瞬態流變學因素,并結合Maxwell粘彈性流變學來模擬板塊震間的粘彈性效應[8]。對于地幔的應力積累和長期演變過程,具有流變學特性的Maxwell粘彈性介質體能很好地模擬不同地塊隨時間緩慢變化的粘彈性流變特性,較好地反映地殼運動過程中產生的地表位移和應力積累過程。同時,Maxwell體也是一種與時間無關的彈性行為和與時間相關的粘性行為的組合,分別以彈性模量和粘度為特征[30]。

本文在模擬2種不同模型的震間板塊運動時采用相同的反向滑動(backslip)設置,即假設一個與同震滑移相反的斷層位錯來模擬不同板塊間相鄰斷裂帶的震間閉鎖[8,31]。大量研究表明,反向滑移模型能較好地解釋震間活動斷裂附近的GPS觀測[8]。本文采用文獻[21]中估計的納斯卡板塊(NAZ)相對于南美板塊(NA)的相對滑動速率作為約束,設置俯沖界面以55 mm/a的相對速率發生滑動,并沿N83°E方向俯沖至南美板塊之下[13,21,32]。

對彈性和粘彈性三維模型,本文設置閉鎖區深度在10～80 km內變化,根據文獻[16]并間隔5 km進行設置。同時為簡化模型,將斷裂閉鎖設置為均勻分布,共完成15個正演模型。使用均方根(root mean square,RMS)準則以量化GPS觀測結果和有限元模型預測值之間的不擬合度,進而通過評估該不擬合度搜索最優的震間閉鎖深度[8]。RMS計算公式為:

(1)

式中,V實測為實測GPS點在東、北方向的速度分量;V模擬為有限元模型在相應點東、北方向的速度分量;n為GPS站點數。

3 結果和討論

圖3為彈性模型與粘彈性模型在不同閉鎖深度的模擬值與GPS實測數據的均方根分布。可以看出,兩種模型表現出相近的變化趨勢:隨著斷層閉鎖深度的增加,彈性模型和粘彈性模型的RMS值先下降,達到最佳擬合后逐漸上升。另外,在斷層閉鎖深度為30～35 km時,彈性模型能較好地擬合GPS數據,RMS值為2.96 mm/a;而對于粘彈性模型,最佳斷層閉鎖深度為25 km,此時對應的RMS值為2.52 mm/a。

圖3 彈性模型與粘彈性模型不同閉鎖深度的模擬值與GPS實測數據的均方根

3.1 最佳彈性模型

采用彈性模型計算的震間最佳閉鎖深度約為30～35 km,與前人研究的秘魯區域閉鎖深度約為30～40 km較為接近[16,21]。由圖4(a)可知,在該閉鎖深度下,彈性模型結果與實測GPS速度場在大部分地區較為接近,但在靠近海溝及北部區域的差異仍較為明顯。從圖4(b)可以看出,除靠近海溝的部分站點外,大部分區域均呈現出NEE向系統性殘差,并與板塊俯沖匯聚方向一致,表明純彈性模型未能較好地模擬研究區震間俯沖變形機制。

圖4 最佳閉鎖深度下GPS數據、模擬值及其殘差

3.2 最佳粘彈性模型

粘彈性模型獲得的最佳閉鎖深度為25 km,與彈性模型相比(RMS為2.96 mm/a),粘彈性模型在整體上能更好地擬合厄瓜多爾地區GPS速度場,模擬值與GPS觀測值的吻合度更高(RMS為2.52 mm/a),結果見圖5。雖然在遠場粘彈性模型的殘差稍大(圖4和5紅色方框區域),但這主要是由粘彈性變形所具有的明顯長波信號所致[8]。此外,遠場區域的構造變形也在一定程度上受安第斯造山運動的影響,同樣也可能導致模擬殘差較大,因此板塊的俯沖匯聚并非其唯一的構造機制[11,16-17]。但從整體上看,相比于純彈性模型,粘彈性模型能明顯提升該區數值模擬結果與GPS數據的擬合度,更好地解釋震間所觀測到的形變場。分析表明,在厄瓜多爾地區,深部粘彈性軟流圈在俯沖板塊的推擠作用下,在震間會引起較為明顯的上地殼變形,且這一變形方向與俯沖匯聚方向基本一致。

圖5 最佳閉鎖深度下GPS數據、模擬值及及其殘差

4 結 語

本文以南美洲板塊高精度GPS速度場為約束,結合區域活動構造特征與深部圈層結構,構建納斯卡板塊與南美洲板塊在厄瓜多爾俯沖區的三維有限元動力學模型,并分別以純彈性與純彈性+粘彈性模型為理論基礎,模擬分析厄瓜多爾地區的震間形變場,估計最佳斷裂閉鎖深度。結果表明,相比于彈性模型,粘彈性模型能更好地解釋震間地表變形,表明深部粘彈性圈層在俯沖板塊的推擠作用下會引起明顯的上地殼變形。此外,基于彈性模型的最佳震間閉鎖深度約為30～35 km,而基于粘彈性模型的最佳閉鎖深度為25 km,進一步說明粘彈性模型同樣會引起俯沖斷裂閉鎖分布特征的差異。今后針對該地區的地球動力學研究有必要進一步考慮深部粘彈性圈層的動力學機制。

本文關注的重點是俯沖界面閉鎖深度對于彈性模型與粘彈性模型在震間影響的差異,并未考慮內陸造山運動與其他斷層構造,在后續研究中會進一步精化模型,研究更復雜情況下的閉鎖深度和區域流變結構。