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

顏復(fù)康 田 鎮(zhèn) 楊志強(qiáng) 楊 兵 梁 沛

1 長安大學(xué)地質(zhì)工程與測繪學(xué)院,西安市雁塔路126號,710054

板塊俯沖是地殼向地幔進(jìn)行物質(zhì)和能量傳輸?shù)闹匾绞街?也是地殼縮短、消亡的主要因素[1-2]。長期的板塊俯沖作用會沿著俯沖界面產(chǎn)生應(yīng)力積累,并在震空區(qū)逐漸增加,當(dāng)應(yīng)力超過巖層承受強(qiáng)度時(shí),則會以地震、海嘯、火山噴發(fā)等形式釋放,給人類生命財(cái)產(chǎn)安全帶來損失[3]。

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

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

1 地質(zhì)背景

南美板塊位于納斯卡板塊、可可斯板塊與大西洋板塊之間,受到來自納斯卡板塊東向強(qiáng)烈沖擊和大西洋海底拉張應(yīng)力,導(dǎo)致板內(nèi)形變非常復(fù)雜[13,18]。

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

2 數(shù)據(jù)和方法

2.1 GPS速度場

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

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

2.2 有限元模型

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

圖2 研究區(qū)有限元模型示意圖

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

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

本文忽略震間模型的重力因素,將其視為絕對應(yīng)力狀態(tài)的擾動(dòng),這些應(yīng)力擾動(dòng)通過沿?cái)鄬咏缑嬉赃\(yùn)動(dòng)學(xué)方式指定反向滑移率進(jìn)行模擬[8]。同時(shí),為使模擬結(jié)果能較客觀地反映該區(qū)域構(gòu)造變形特征,建模時(shí)忽略地表曲率[27]。有限元節(jié)點(diǎn)上的累積位移可根據(jù)其在假定笛卡爾坐標(biāo)系內(nèi)的原始位置進(jìn)行計(jì)算。即使在遠(yuǎn)場中,地球曲率導(dǎo)致的水平差異也可以忽略不計(jì),球面坐標(biāo)系和笛卡爾坐標(biāo)系測試模型之間的差異很小[23,27]。此外,已有研究表明,在模擬百年尺度的震間變形時(shí),該區(qū)域歷史大地震的影響較小,因此本文在之后研究中并未考慮歷史地震的影響[8]。

2.3 巖石圈流變結(jié)構(gòu)配置

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

表1 地塊屬性

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

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

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

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

(1)

式中,V實(shí)測為實(shí)測GPS點(diǎn)在東、北方向的速度分量;V模擬為有限元模型在相應(yīng)點(diǎn)東、北方向的速度分量;n為GPS站點(diǎn)數(shù)。

3 結(jié)果和討論

圖3為彈性模型與粘彈性模型在不同閉鎖深度的模擬值與GPS實(shí)測數(shù)據(jù)的均方根分布?？梢钥闯?兩種模型表現(xiàn)出相近的變化趨勢:隨著斷層閉鎖深度的增加,彈性模型和粘彈性模型的RMS值先下降,達(dá)到最佳擬合后逐漸上升。另外,在斷層閉鎖深度為30～35 km時(shí),彈性模型能較好地?cái)M合GPS數(shù)據(jù),RMS值為2.96 mm/a;而對于粘彈性模型,最佳斷層閉鎖深度為25 km,此時(shí)對應(yīng)的RMS值為2.52 mm/a。

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

3.1 最佳彈性模型

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

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

3.2 最佳粘彈性模型

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

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

4 結(jié) 語

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

本文關(guān)注的重點(diǎn)是俯沖界面閉鎖深度對于彈性模型與粘彈性模型在震間影響的差異,并未考慮內(nèi)陸造山運(yùn)動(dòng)與其他斷層構(gòu)造,在后續(xù)研究中會進(jìn)一步精化模型,研究更復(fù)雜情況下的閉鎖深度和區(qū)域流變結(jié)構(gòu)。