王迎 李江海,? 馬昌明 宋玨琛

基于離散元法的龍門(mén)山構(gòu)造帶隆升變形主控因素研究

王迎1,2李江海1,2,?馬昌明1,2宋玨琛1,2

1.造山帶與地殼演化教育部重點(diǎn)實(shí)驗(yàn)室, 北京大學(xué)地球與空間科學(xué)學(xué)院, 北京 100871; 2.北京大學(xué)石油與天然氣研究中心, 北京 100871; ?通信作者, E-mail: jhli@pku.edu.cn

為探明龍門(mén)山構(gòu)造帶隆升變形的主要控制因素, 基于龍門(mén)山構(gòu)造帶東西兩側(cè)下地殼物質(zhì)層屬性差異巨大的特征, 進(jìn)行 3 組 PFC2D 離散元數(shù)值模擬對(duì)比實(shí)驗(yàn), 將深度擴(kuò)大至下地殼, 記錄顆粒運(yùn)動(dòng)狀態(tài), 實(shí)現(xiàn)定量化分析。實(shí)驗(yàn)得到的變形結(jié)果及模型顆粒運(yùn)動(dòng)矢量圖顯示, 在下地殼物質(zhì)屬性無(wú)明顯差異的條件下, 板塊碰撞擠壓應(yīng)力及地殼厚度的差異不會(huì)在龍門(mén)山構(gòu)造帶形成巨大的地形高差。當(dāng)下地殼黏度系數(shù)存在明顯差異時(shí), 軟弱下地殼物質(zhì)層顆粒相對(duì)運(yùn)動(dòng)速率為 1.5～2.94m/s, 平均運(yùn)動(dòng)速率為 1.62m/s, 大約是堅(jiān)硬下地殼層顆粒平均運(yùn)動(dòng)速率的 54 倍。模型中部(龍門(mén)山構(gòu)造帶)出現(xiàn)隆升變形, 縱向影響范圍為 94.74%, 隆升幅度為19.85%。軟弱下地殼上覆的中地殼和上地殼顆粒具有較大的向上速度分量, 上地殼物質(zhì)層上涌趨勢(shì)明顯。巴顏喀拉塊體和四川盆地地殼存在 20km 的厚度差異, 使得龍門(mén)山構(gòu)造帶隆升幅度由 14.79% 增至 19.85%。綜合分析 3 組離散元模擬實(shí)驗(yàn)結(jié)果, 得出巴顏喀拉地塊下地殼物質(zhì)層與四川盆地下地塊物質(zhì)層的黏度差異是龍門(mén)山構(gòu)造帶垂向隆升變形最關(guān)鍵控制因素的結(jié)論, 在下地殼黏度結(jié)構(gòu)存在明顯差異的前提下, 巴顏喀拉塊體和四川盆地的地殼厚度差異對(duì)龍門(mén)山構(gòu)造帶縱向上逆沖隆升幅度有明顯的促進(jìn)作用。

龍門(mén)山構(gòu)造帶; 離散元數(shù)值模擬; 下地殼異質(zhì)性; 地殼厚度; 主控因素

龍門(mén)山逆沖構(gòu)造帶是青藏高原東部巴顏喀拉地塊與四川盆地的地質(zhì)、地貌和地球物理邊界, 涉及印支造山期中國(guó)大陸主體拼合以及喜馬拉雅造山期青藏高原隆升兩大構(gòu)造事件[1?3]。在經(jīng)歷破壞力極大的 2008 年汶川地震后, 作為最活躍的構(gòu)造邊界和青藏高原周緣生長(zhǎng)的關(guān)鍵區(qū)域, 龍門(mén)山逆沖構(gòu)造帶成為全球構(gòu)造地質(zhì)學(xué)家研究逆沖構(gòu)造隆升機(jī)制的熱點(diǎn)地區(qū)。關(guān)于龍門(mén)山逆沖構(gòu)造帶隆升變形機(jī)制的討論集中在地殼縮短(crustal shortening)[4?5]和下地殼管道流(poiseuille flow)[6?8]兩大機(jī)制。地殼縮短機(jī)制認(rèn)為, 板塊碰撞擠壓應(yīng)力引起上地殼在水平方向上發(fā)生縮短, 由此產(chǎn)生的斷層及其相關(guān)褶皺構(gòu)造和逆沖推覆構(gòu)造是龍門(mén)山地區(qū)地殼加厚隆升和地形高差巨大的主要原因。下地殼管道流機(jī)制認(rèn)為, 青藏高原深部存在黏度低的易發(fā)生韌性流變變形的物質(zhì)層, 其東南向流動(dòng)受到四川盆地阻擋, 在龍門(mén)山地區(qū)形成陡峭的地形。前者將研究重點(diǎn)集中在塊體間的相互作用方面, 缺少對(duì)塊體內(nèi)部運(yùn)動(dòng)和應(yīng)力狀態(tài)的分析; 后者則是在理想狀態(tài)下建立流變學(xué)模型, 忽略地殼縮短量以及青藏高原隆升所產(chǎn)生的高地形重力載荷的影響。

前人通過(guò)低溫?zé)崮甏鷮W(xué)、大地電磁探測(cè)和地震動(dòng)力學(xué)數(shù)值模擬等手段, 得出地殼縮短[9]、地殼黏度結(jié)構(gòu)[10?12]和地殼厚度差異[12?13]是晚新生代龍門(mén)山構(gòu)造隆升變形重要影響因素的結(jié)論。本文基于前人對(duì)該地區(qū)的地質(zhì)和地球物理研究成果, 利用離散元模擬技術(shù), 對(duì)顆粒的運(yùn)動(dòng)狀態(tài)進(jìn)行記錄和定量分析, 獲取更多、更準(zhǔn)確的動(dòng)力學(xué)信息, 通過(guò)單因素變量控制方法, 探討擠壓背景下的下地殼黏度差異和地殼厚度對(duì)龍門(mén)山構(gòu)造帶隆升變形的影響, 分析龍門(mén)山構(gòu)造帶的隆升變形機(jī)理。

1 區(qū)域地質(zhì)背景

龍門(mén)山構(gòu)造帶位于特提斯?喜馬拉雅構(gòu)造域與濱西太平洋構(gòu)造域的交接轉(zhuǎn)換部位[14?15], 東側(cè)以安縣?灌縣斷裂與四川盆地分界, 西側(cè)以茂縣?汶川斷裂與巴顏喀拉地塊分界, 南部為青藏高原東南緣的川滇構(gòu)造帶, 北部為秦嶺造山帶南緣的米倉(cāng)山構(gòu)造帶[16](圖 1(a))。龍門(mén)山構(gòu)造帶是典型的陸內(nèi)擠壓造山帶, 晚三疊世之前處于被動(dòng)大陸邊緣拉張背景下, 晚三疊世末期轉(zhuǎn)變?yōu)闃?gòu)造擠壓環(huán)境。受印支期造山運(yùn)動(dòng)以及喜馬拉雅造山運(yùn)動(dòng)的影響, 巴顏喀拉塊體變形劇烈, 龍門(mén)山構(gòu)造帶發(fā)生強(qiáng)烈的沖斷隆升[17]。龍門(mén)山褶皺沖斷帶平均海拔高度為 2000～4000km (圖 1(b)), 其主峰九頂山海拔高度達(dá)到 4989km。沿NE-SW 走向, 龍門(mén)山構(gòu)造帶北段地形梯度小于中段和南段[18]。在研究區(qū)范圍內(nèi), 龍門(mén)山構(gòu)造斷裂帶位于莫霍面深度變化最為強(qiáng)烈的地帶, 自西北向東南方向, 巴顏喀拉地塊和四川盆地之間存在巨大的莫霍面深度差(圖 1(c))。

龍門(mén)山構(gòu)造帶的總體長(zhǎng)度約為 600km, 自西向東發(fā)育 3 條顯露于地表的主干沖斷裂: 青川?茂汶斷裂、北川?映秀斷裂和安縣?灌縣斷裂[19?20]。青川?茂汶斷裂為龍門(mén)山構(gòu)造帶與松潘?甘孜褶皺帶的邊界斷裂[21]; 北川?映秀斷裂為汶川大地震的發(fā)震斷裂, 其地表破裂帶長(zhǎng)度為 240～300km[22]; 安縣?灌縣斷裂為其南東邊界。晚新生代以來(lái), 它們都有明顯的活動(dòng)跡象, 且具有發(fā)生強(qiáng)震的可能性。它們?cè)诖瓜蛏铣石B瓦狀展布, 向四川盆地方向逆沖推覆, 在地下 20km 深處匯聚成一條大型剪切帶[23?24]。

(a)研究區(qū)地形及斷裂分布; (b)龍門(mén)山斷裂帶地形剖面; (c)龍門(mén)山斷裂帶莫霍面深度。地形數(shù)據(jù)來(lái)源于地理空間數(shù)據(jù)云(http://www.gscloud. cn), 地震數(shù)據(jù)來(lái)源于中國(guó)地震臺(tái)網(wǎng)(https://news.ceic.ac.cn), 莫霍面深度數(shù)據(jù)來(lái)源于Crust1.0 (https://igppweb.ucsd.edu/～gabi/crust1.html)

2 離散元數(shù)值模擬

2.1 離散元方法及實(shí)驗(yàn)參數(shù)

Cundall 等[25]于 1979 年提出離散元(discrete ele-ment merhod)方法的理論基礎(chǔ)。該數(shù)值模擬方法基于顆粒單元間的接觸準(zhǔn)則, 利用時(shí)間步長(zhǎng)和有限差分方法求解系統(tǒng)內(nèi)每個(gè)顆粒單元的牛頓運(yùn)動(dòng)學(xué)方程[25?26]。求解時(shí), 將地質(zhì)目標(biāo)體離散為質(zhì)量有限的圓盤(pán)狀的剛性顆粒, 顆粒間為柔性接觸, 接觸范圍極小(點(diǎn)接觸), 接觸位置可發(fā)生重疊, 且存在特殊的鏈接強(qiáng)度。通過(guò)調(diào)整顆粒強(qiáng)度、顆粒之間的黏結(jié)強(qiáng)度以及內(nèi)摩擦系數(shù)等微觀參數(shù), 設(shè)定和模擬不同的巖石類(lèi)型和巖性組合?；趧?dòng)態(tài)應(yīng)力松弛顯式差分法求解運(yùn)動(dòng)方程和動(dòng)力方程, 采用搜索算法搜索顆粒體接觸條件, 并計(jì)算接觸受力狀態(tài), 將顆粒單元間的不平衡力重新作用到節(jié)點(diǎn)上, 迭代至整個(gè)系統(tǒng)不平衡力足夠小或者節(jié)點(diǎn)位移趨于平衡為止, 達(dá)到有效模擬自然重力沉積環(huán)境下非連續(xù)介質(zhì)大應(yīng)變量構(gòu)造變形的目的[27?28]。近年來(lái), 離散元數(shù)值模擬方法廣泛應(yīng)用于構(gòu)造地質(zhì)研究, 如伸展變形[29]、板塊俯沖變形[30]和褶皺沖斷帶構(gòu)造變形[31]等。本文利用 PFC (particle flow code) 2D 軟件進(jìn)行離散元數(shù)值模擬。

2.1.1幾何參數(shù)

隨著觀測(cè)數(shù)據(jù)的增加和地球物理學(xué)方法(接收函數(shù)、面波層析成像、體波層析成像和深地震測(cè)深等)的進(jìn)步, 研究區(qū)地殼結(jié)構(gòu)模型愈發(fā)精細(xì)。研究表明, 龍門(mén)山地處莫霍面向西傾斜的陡變帶, 其西側(cè)巴顏喀拉塊體的上地殼、中地殼和下地殼厚度分別為 20～22, 26～30 和 14～18km, 地殼總厚度達(dá)到60～70km; 其東側(cè)四川盆地的上地殼、中地殼和下地殼厚度分別為 10～12, 10～15 和 15～20km, 地殼總厚度僅為 30～40km[32?34]。依據(jù)研究區(qū)實(shí)際地殼結(jié)構(gòu), 本文離散元數(shù)值模型與實(shí)際地質(zhì)體的尺寸按 1:1000 的比例設(shè)置, 長(zhǎng) 400m, 寬 100m。為模擬莫霍面的西傾形態(tài), 將模型底部設(shè)置為非水平邊界, 在軸的 100～275m 區(qū)間為傾斜邊界, 傾角大約為6°。墻體法向剛度和切向剛度均設(shè)置為 1010N/m。在 400m′100m 區(qū)域內(nèi)生成隨機(jī)粒子, 顆粒半徑在0.8～1.0m 之間, 不同粒徑的顆粒數(shù)量服從高斯分布。加載重力加速度為 9.81m/s2, 允許顆粒聚集, 自然壓實(shí)。達(dá)到平衡狀態(tài)時(shí), 刪除頂部多余顆粒, 僅保留 400m′60m范圍內(nèi)顆粒。

2.1.2微觀物理學(xué)參數(shù)

本文研究區(qū)下地殼實(shí)際黏度系數(shù)高達(dá) 1019～1021Pa·s[35?37], 如果在離散元模型中如此設(shè)定, 單個(gè)模型的模擬運(yùn)行時(shí)間會(huì)長(zhǎng)達(dá)數(shù)年。本文參照前人的模型設(shè)定方法以及宏觀參數(shù)與微觀參數(shù)關(guān)系測(cè)定結(jié)果[38?39], 并進(jìn)行大量參數(shù)標(biāo)定, 在滿(mǎn)足模擬需求和保證模擬有效性的前提下, 為了提高計(jì)算效率, 將二維離散模型中下地殼物質(zhì)層顆粒的黏度系數(shù)縮小1010倍, 設(shè)定為 109～1011Pa·s。

合理的參數(shù)是保證數(shù)值模擬實(shí)驗(yàn)結(jié)果有效性的前提。在通過(guò)二維顆粒流分析程序(PFC2D)進(jìn)行離散元數(shù)值模擬實(shí)驗(yàn)時(shí), 多通過(guò)雙軸力學(xué)模擬實(shí)驗(yàn)的方法來(lái)確定顆粒的微觀參數(shù)[40?41]。從前人的研究結(jié)果中獲得的研究區(qū)巖石力學(xué)參數(shù)如下: 巴顏喀拉塊體以四川盆地上地殼物質(zhì)層的內(nèi)聚力分別為 30 和50MPa[42], 內(nèi)摩擦角分別為 10°和 30°[43], 楊氏模量均為 8.75×1010Pa[44], 泊松比均為 0.25[45]。巴顏喀拉地塊以及四川盆地中下地殼物質(zhì)層的楊氏模量均為1.0×1011Pa[44], 泊松比均為 0.25[45]。以上述巖石力學(xué)參數(shù)為基礎(chǔ), 進(jìn)行雙軸力學(xué)數(shù)值模擬實(shí)驗(yàn)調(diào)試, 施加軸向恒定壓縮速率, 平臺(tái)伺服施加指定圍壓, 利用 Servo 系統(tǒng)和 Get_gain 函數(shù)獲得巖石試樣破壞全應(yīng)力?應(yīng)變曲線(xiàn)。通過(guò)圍壓?最大軸主應(yīng)力對(duì)應(yīng)的應(yīng)力摩爾圓及包絡(luò)線(xiàn), 對(duì)比實(shí)際的材料參數(shù), 確定模型微觀輸入?yún)?shù), 獲得各層 PFC2D 微觀參數(shù)。本文采用雙軸數(shù)值試驗(yàn)設(shè)定實(shí)驗(yàn)參數(shù), 并參考文獻(xiàn)[30,40,46? 49]。模型各層微觀物理學(xué)參數(shù)見(jiàn)表 1。

2.1.3邊界條件

地表 GPS 速度觀測(cè)數(shù)據(jù)顯示, 印度板塊以 40mm/a 的速度沿喜馬拉雅山脈向歐亞板塊推進(jìn), 青藏高原受到強(qiáng)烈擠壓發(fā)生連續(xù)變形, 中部的巴顏喀拉塊體整體呈東南向運(yùn)動(dòng), 速率約為 21mm/a[50?51]。依據(jù)動(dòng)力學(xué)相似原則設(shè)定實(shí)驗(yàn)?zāi)Ｐ妥冃嗡俾? 則Ramberg (m)系數(shù)(無(wú)量綱)的計(jì)算公式[52]如下:

式中,,,,和分別表示密度、重力加速度、幾何長(zhǎng)度、黏度以及變形速率, 下角標(biāo)m和n分別代表實(shí)驗(yàn)?zāi)Ｐ秃妥匀唤绲刭|(zhì)原型。在理想條件下, 地質(zhì)原型和實(shí)驗(yàn)?zāi)Ｐ偷膍值應(yīng)相同或?qū)儆谕涣考?jí)。將本文地質(zhì)原型及表 2 中參數(shù)代入式(1), 得到實(shí)驗(yàn)?zāi)Ｐ妥髠?cè)的擠壓速率為 6×10?3mm/s。

2.2 實(shí)驗(yàn)設(shè)計(jì)

2.2.1黏度結(jié)構(gòu)對(duì)地表變形的影響

在龍門(mén)山構(gòu)造帶深部, 黏度呈現(xiàn)強(qiáng)烈的橫向變化, 東西兩側(cè)下地殼存在 1～2 個(gè)量級(jí)的黏度差異。巴顏喀拉塊體下地殼黏度范圍為 1.48′1017～1.12′1021Pa·s, 平均值為 1.0′1019Pa·s; 四川盆地下地殼黏度范圍為 4.09×1019～7.08×1020Pa·s, 平均值為 1.0×1020Pa·s[35?37]。為考察下地殼黏度結(jié)構(gòu)對(duì)龍門(mén)山構(gòu)造帶隆升變形的影響, 本文設(shè)計(jì)兩組對(duì)比模型: 模型 1和模型 2。

模型 1 將 400m′60m 范圍內(nèi)的隨機(jī)粒子劃分為 4 個(gè)顆粒集, 自上而下分別為上地殼、中地殼和下地殼, 各層顆粒沉積密度均沿軸正方向減小。下地殼分為左右兩段, 左段軟弱下地殼的顆粒黏度系數(shù)設(shè)為 109Pa·s, 右段堅(jiān)硬下地殼的顆粒黏度系數(shù)設(shè)為 1010Pa·s, 分別對(duì)應(yīng)巴顏喀拉地塊下地殼和四川盆地下地殼(圖 2(a))。上地殼和中地殼內(nèi)設(shè)標(biāo)志層, 便于后期識(shí)別構(gòu)造變形。為了保證力和力矩的有效傳遞, 下地殼層顆粒接觸類(lèi)型為線(xiàn)性模型, 上地殼層和中地殼層顆粒接觸類(lèi)型為線(xiàn)性平行黏結(jié)模型。模型右側(cè)固定, 左側(cè)以 6′10?3mm/s的速度擠壓。記錄擠壓推進(jìn) 2.5%, 5%和 7.5%時(shí)的變形情況(圖 3)。

表1 離散元模擬微觀物理學(xué)參數(shù)

表2 實(shí)驗(yàn)材料參數(shù)及模型相似比

說(shuō)明: 模型各層密度與自然界地質(zhì)體實(shí)際密度一致; 黏度數(shù)值參考文獻(xiàn)[35,53?54]。

模型 2 僅在縱向上做顆粒集的分層劃分, 將400m′60m 范圍內(nèi)的隨機(jī)粒子分為 3 個(gè)顆粒集(圖2(b))。下地殼物質(zhì)層黏度無(wú)差別, 為 1010Pa·s。其他設(shè)置(地殼厚度、推進(jìn)速率和顆粒接觸類(lèi)型)均與模型 1 相同。記錄擠壓推進(jìn) 2.5%, 5%和 7.5%時(shí)的變形情況(圖 4)。

2.2.2地殼厚度差異對(duì)地表變形的影響

在印度洋板塊的北向擠壓下, 青藏高原具有明顯的 NNE 向運(yùn)動(dòng)趨勢(shì), 受西伯利亞板塊和塔里木地塊的阻擋, 發(fā)生快速隆起, 其中部的巴顏喀拉塊體成為巨厚陸殼體[55?56]。

為探討巴顏喀拉塊體與四川盆地地殼厚度差異對(duì)龍門(mén)山構(gòu)造帶隆升變形的影響, 本文設(shè)計(jì)模型 3, 與模型 1 進(jìn)行對(duì)比。

模型 3 保留 400m′40m范圍內(nèi)的顆粒, 分為 4 個(gè)顆粒集。自上而下劃分為上地殼、中地殼和下地殼 3 層。下地殼劃分為軟弱下地殼和堅(jiān)硬下地殼, 黏度系數(shù)均與模型 1 一致, 各物質(zhì)層顆粒沉積密度沿軸正方向保持不變。模型 3中地殼厚度均一, 即莫霍面不存在起伏, 因此模型3 的底板設(shè)置為水平邊界(圖 2(c))。其他設(shè)置(黏度系數(shù)、推進(jìn)速率和顆粒接觸類(lèi)型)均與模型 1 相同。記錄擠壓推進(jìn) 2.5%, 5.0%和 7.5%時(shí)的變形情況(圖 5)。

3 離散元數(shù)值模擬結(jié)果

提取二維離散元模擬變形結(jié)果(圖 3～5)及顆粒運(yùn)動(dòng)矢量圖(圖 6), 后者記錄了模型顆粒在末位置處的累積速度。模型 2 和 3 在近擠壓端出現(xiàn)速度場(chǎng)奇點(diǎn)異?，F(xiàn)象(圖 6(b)和(c)中紅色箭頭所示), 奇點(diǎn)數(shù)量占比低于 0.05%, 對(duì)模型整體變形趨勢(shì)及速度場(chǎng)的影響可忽略不計(jì)。

如圖 3 所示, 模型 1 中, 隨著左端的持續(xù)擠壓, 當(dāng)應(yīng)力超過(guò)一定范圍時(shí), 顆粒在互斥作用下發(fā)生剪切錯(cuò)動(dòng)和分離, 宏觀上表現(xiàn)為地層垂向的逆沖抬升和水平方向的縮短變形。在擠壓量達(dá)到 2.5%的早期擠壓變形階段, 左段黏度小的軟弱下地殼受左側(cè)擠壓作用發(fā)生橫向流動(dòng), 帶動(dòng)其上部物質(zhì)層同步右移, 受到右段黏度大的不易流動(dòng)的堅(jiān)硬下地殼阻擋, 出現(xiàn)堆積聚集現(xiàn)象, 在中部(下地殼物質(zhì)層屬性變化處)出現(xiàn)褶曲, 高角度逆沖斷裂 F1和 F2相繼形成。F2與隨后形成的反向調(diào)節(jié)斷層 F3構(gòu)成“V”字形的斷層組合, 模型左部形成正斷層 F4。在擠壓量達(dá)到5.0%的中期擠壓變形階段, 模型左段斷層數(shù)量增多, 靠近擠壓端形成沖斷斷層 F7, 模型左段和右段差異變形態(tài)勢(shì)逐漸增強(qiáng)。擠壓量達(dá)到 7.5%時(shí), 斷層數(shù)量未增加, F1, F2和F6的斷距小幅度增加。

圖2 離散元模擬實(shí)驗(yàn)的初始模型設(shè)計(jì)

圖3 離散元模擬實(shí)驗(yàn)?zāi)Ｐ?變形過(guò)程

圖4 離散元模擬實(shí)驗(yàn)?zāi)Ｐ?變形過(guò)程

圖5 離散元模擬實(shí)驗(yàn)?zāi)Ｐ?變形過(guò)程

圖6 離散元模擬實(shí)驗(yàn)顆粒運(yùn)動(dòng)矢量對(duì)比

如圖 4 所示, 模型 2 中右段發(fā)育褶曲疊瓦狀逆沖斷裂, 而在物質(zhì)層厚度相對(duì)較大的左段和中部未發(fā)生明顯的變形。模型顆粒運(yùn)動(dòng)速率整體上較小, 為 0～0.005m/s。沿軸正方向, 模型顆粒的向上速度分量逐漸增大, 遠(yuǎn)擠壓端顆粒顯示較強(qiáng)的上涌趨勢(shì)。

如圖 5 所示, 模型 3 中的變形集中在黏度系數(shù)較小的軟弱下地殼區(qū)段, 近擠壓端出現(xiàn)反沖斷層, 模型中部(下地殼物質(zhì)層屬性變化處)形成變形強(qiáng)烈的隆起構(gòu)造。軟弱下地殼層顆粒具有較大的水平速度分量, 相對(duì)運(yùn)動(dòng)速率為 1.0～1.75m/s, 其上覆中地殼和上地殼物質(zhì)層顆粒具有較大的垂向速度分量, 相對(duì)運(yùn)動(dòng)速率為 1.25～2.25m/s。堅(jiān)硬下地殼上覆中地殼和上地殼物質(zhì)層顆粒運(yùn)動(dòng)速率趨近零。

4 討論

模型 1 左右兩段下地殼的黏度和厚度存在明顯差異, 左側(cè)推進(jìn)擠壓作用導(dǎo)致模型中部(龍門(mén)山構(gòu)造帶)發(fā)生隆升變形, 縱向影響范圍為 94.74%, 隆升幅度為 19.85%。以下地殼物質(zhì)層屬性變化處為界, 其左右兩側(cè)的模型顆粒運(yùn)動(dòng)方向和相對(duì)運(yùn)動(dòng)速率存在較大的差異, 右側(cè)顆粒僅在有限空間范圍內(nèi)做低速運(yùn)動(dòng), 左側(cè)顆粒的相對(duì)運(yùn)動(dòng)速率整體上較大, 且在水平方向上具有明顯的分段性。軟弱下地殼物質(zhì)層顆粒具有較大的水平運(yùn)動(dòng)速度分量, 上地殼和中地殼物質(zhì)層具有較大的縱向運(yùn)動(dòng)速度分量(圖 3和圖 6(a))。

上述結(jié)果與前人通過(guò)分析該地區(qū)的重力異常數(shù)據(jù)、三維成像反演結(jié)果及 GPS 水平位移速度等資料得出的認(rèn)識(shí)匹配度良好。在歐亞板塊框架下, 巴顏喀拉地塊的構(gòu)造運(yùn)動(dòng)特征較為復(fù)雜, 主要處于NW-SE向拉張和NE-SW向擠壓變形狀態(tài), 塊體內(nèi)部的運(yùn)動(dòng)速率及其差異性明顯, 自西向東運(yùn)動(dòng)速率逐漸減小[57]。中國(guó)地殼運(yùn)動(dòng)觀測(cè)網(wǎng)絡(luò)數(shù)據(jù)[58]顯示, 巴顏喀拉塊體中部運(yùn)動(dòng)速率約為 20mm/a, 東部約為 10mm/a。巴顏喀拉地塊是青藏高原深部軟弱物質(zhì)發(fā)生東南向“逃逸”及塊體相互作用、變形的主要承載體[57?69]。四川盆地深部存在高速度、低電導(dǎo)率的下地殼物質(zhì)層, 對(duì)青藏高原軟弱下地殼的東南向“逃逸”起攔截作用[4,60]。

模型 2 中地殼厚度存在明顯差異, 左右兩段下地殼的黏度系數(shù)相同。在模型左側(cè)的擠壓下, 變形出現(xiàn)在物質(zhì)層厚度較薄的右端, 全空間區(qū)域顆粒相對(duì)運(yùn)動(dòng)速率均小于 1m/s。綜合分析模型 1 和 2 的變形結(jié)果(圖 3 和 4)和顆粒運(yùn)動(dòng)速度場(chǎng)(圖 5)可知, 巴顏喀拉地塊下地殼物質(zhì)層與四川盆地下地殼物質(zhì)層黏度結(jié)構(gòu)的差異是龍門(mén)山構(gòu)造帶縱向隆升變形的先決條件。在下地殼黏度結(jié)構(gòu)沒(méi)有明顯差異的前提下, 來(lái)自巴顏喀拉塊體的南東向擠壓應(yīng)力以及龍門(mén)山構(gòu)造帶東西兩側(cè)地殼厚度的差異, 對(duì)龍門(mén)山構(gòu)造帶縱向地貌形態(tài)的塑造作用不明顯。

模型 3 中, 左右兩段下地殼的黏度存在明顯的差異, 地殼厚度相同。在左側(cè)的擠壓作用下, 近擠壓端發(fā)育反沖斷層, 模型中部(龍門(mén)山構(gòu)造帶)發(fā)育構(gòu)造三角帶, 縱向影響范圍為 93.17%, 隆升幅度為14.79%。

綜合分析模型 1 和模型 3 的變形結(jié)果(圖 3 和 5)及顆粒運(yùn)動(dòng)速度場(chǎng)(圖 6)可知, 巴顏喀拉塊體與四川盆地地殼厚度的差異對(duì)龍門(mén)山構(gòu)造帶垂向隆升變形起到一定的促進(jìn)作用。在下地殼黏度結(jié)構(gòu)存在明顯差異的前提下, 龍門(mén)山東西兩側(cè)地殼厚度相同時(shí), 隆升幅度僅為 14.79%; 龍門(mén)山東西兩側(cè)地殼存在20km 的厚度差時(shí), 隆升幅度可以達(dá)到 19.85%。與模型 3 相比, 模型 1 左側(cè)地殼、中地殼以及下地殼物質(zhì)層顆粒相對(duì)運(yùn)動(dòng)速率在水平方向具有明顯的分段性, 縱向上具有較高的一致性。由此推測(cè), 地殼厚度差異產(chǎn)生的壓力梯度能夠促進(jìn)軟弱下地殼的水平運(yùn)動(dòng), 進(jìn)而影響上、中地殼的運(yùn)動(dòng)狀態(tài)。在下地殼物質(zhì)層屬性變化處, 模型 1 和 3 的顆粒相對(duì)運(yùn)動(dòng)速率及變形趨勢(shì)均急劇降低, 說(shuō)明龍門(mén)山東西兩側(cè)下地殼黏度的差異是其構(gòu)造帶縱向隆升變形的先決條件。

5 結(jié)論

本文基于 PFC2D 離散元數(shù)值模擬對(duì)比實(shí)驗(yàn), 針對(duì)地殼黏度結(jié)構(gòu)和地殼厚度, 研究龍門(mén)山構(gòu)造帶隆升變形的影響作用, 得到以下主要結(jié)論。

1)下地殼異質(zhì)性是龍門(mén)山構(gòu)造帶縱向隆升變形的最為重要的控制因素。在巴顏喀拉地塊與四川盆地下地殼物質(zhì)層黏度系數(shù)沒(méi)有明顯差異的情況下, 地殼厚度的差異不會(huì)形成龍門(mén)山的縱向高地貌形態(tài)。

2)在下地殼黏度結(jié)構(gòu)存在明顯差異的前提下, 巴顏喀拉塊體與四川盆地的地殼厚度的差異對(duì)龍門(mén)山構(gòu)造帶的縱向逆沖隆升幅度有明顯的促進(jìn)作用。龍門(mén)山東西兩側(cè)地殼厚度存在 20km 的厚度差時(shí), 隆升幅度可以達(dá)到 19.85%; 龍門(mén)山東西兩側(cè)地殼厚度相同時(shí), 隆升幅度僅為 14.79%。

[1]Clark M K, Bush J W M, Royden L H. Dynamic topography produced by lower crustal flow. against rheological strength heterogeneities bordering the Ti-betan Plateau. Geophysical Journal International, 2005, 162(2): 575?590

[2]Harrowfield M J, Wilson C J L. Indosinian deforema-tion of the Songpan-Garze Fold Belt. northeast Tibe- tan Plateau. Journal of Structural Geology, 2005, 27: 101?117

[3]Wilson C J L, Harrowfield M J, Reid A J. Brittle modification of Triassic architecture in eastern. Tibet: implications for the construction of the Cenozoic plateau. Journal of Asian Earth Sciences, 2006, 27(3): 341?357

[4]Hubbard J, Shaw J H. Uplift of the Longmen Shan and Tibetan plateau, and the 2008 Wenchuan (=7.9) earthquake. Translated World Ssmology, 2009, 458: 194?197

[5]Tapponnier P, Xu Z Q, Roger F, et al. Oblique step-wise rise and growth of the Tibet Plateau. Science, 2001, 294: 1671?1677

[6]Clark M K, Royden L H. Topographic ooze: building the eastern margin of Tibet by lower crustal flow. Geology, 2000, 28(8): 703?706

[7]Cook K L, Royden L H. The role of crustal strength variations in shaping collisional orogens, with app-lication to Tibet. Journal of Geophysical Research: Solid Earth, 2008, 113(B8): B08407

[8]顏丹平, 孫銘, 鞏凌霄, 等. 青藏高原東緣龍門(mén)山前陸逆沖帶復(fù)合結(jié)構(gòu)與生長(zhǎng). 地質(zhì)力學(xué)學(xué)報(bào), 2020, 26(5): 615?633

[9]譚錫斌, 徐錫偉, 李元希, 等. 龍門(mén)山中段中央斷裂和前山斷裂的晚新生代垂向活動(dòng)性差異及其構(gòu)造意義. 地球物理學(xué)報(bào), 2015, 58(1): 143?152

[10]朱濤, 詹艷, 孫翔宇, 等. 四川龍門(mén)山斷裂帶高精度地殼/巖石圈黏度結(jié)構(gòu)及其動(dòng)力學(xué)意義. 地球物理學(xué)報(bào), 2020, 63(1): 196?209

[11]Liu C, Zhu B J, Yang X L, et al. Crustal rheology control on earthquake activity across the eastern mar-gin of the Tibetan Plateau: insights from numerical modelling. Journal of Asian Earth Sciences, 2015, 100: 20?30

[12]England P, Houseman G. Finite strain calculations of continental defornation. 2. comparison with the India-Asia collision zone. Journal of Geophysical Research: Solid Earth and Planets, 1986, 91(B3): 3664?3676

[13]李廷棟. 青藏高原隆升的過(guò)程和機(jī)制. 地球?qū)W報(bào), 1995(1): 1?9

[14]Enkin R J, Yang Z Y, Chen Y, et al. Paleomagnetic constraints on the geodynamic: history of the Major Blocks of China from the permian to the present. Journal of Geophysical Research: Solid Earth, 1992, 97(B10): 13953?13989

[15]Yin A, Harrison T M. Geologic evolution of the Himalayan-Tibetan Orogen. Annual Review of Earth and Planetary Sciences, 2000, 28: 211?280

[16]李智武, 宋天慧, 王自劍, 等. 川西?龍門(mén)山盆山系統(tǒng)走向差異演化的變形、隆升和沉積記錄及關(guān)鍵構(gòu)造變革期討論. 成都理工大學(xué)學(xué)報(bào)(自然科學(xué)版), 2021, 48(3): 257?282

[17]Jia D, Wei G, Chen Z, et al. Longmen Shan fold-thrust belt and its relation to the western Sichuan Basin in central China: new insights from hydrocar-bon exploration. AAPG Bulletin, 2006, 90(9): 1425? 1447

[18]鄧賓, 何宇, 黃家強(qiáng), 等. 龍門(mén)山褶皺沖斷帶擴(kuò)展生長(zhǎng)過(guò)程——基于低溫?zé)崮甏鷮W(xué)模型證據(jù). 地質(zhì)學(xué)報(bào), 2019, 93(7): 1588?1600

[19]金文正, 湯良杰, 楊克明, 等. 川西龍門(mén)山褶皺 沖斷帶分帶性變形特征. 地質(zhì)學(xué)報(bào), 2007, 81(8): 1072?1080

[20]Jin W, Tang L, Yang K, et al. Segmentation of the Longmen Mountains thrust belt, Western Sichuan Foreland Basin, SW China. Tectonophysics, 2010, 485: 107?121

[21]鄭勇, 李海兵, 王煥, 等. 印支期龍門(mén)山斷裂帶的逆沖?推覆構(gòu)造和沉積響應(yīng). 地質(zhì)論評(píng), 2018, 64(1): 45?61

[22]李海兵, 王宗秀, 付小方, 等. 2008 年 5 月 12 日汶川地震(Ms 8.0)地表破裂帶的分布特征. 中國(guó)地質(zhì), 2008, 35(5): 803?813

[23]何祥麗, 李海兵, 王煥, 等. 蠕滑斷裂帶巖石組成和構(gòu)造特征分析: 以龍門(mén)山灌縣?安縣斷裂帶為例. 巖石學(xué)報(bào), 2020, 36(10): 299?314

[24]Burchfiel B C, Royden L H, Van der Hilst R D, et al. A geological and geophysical context for the Wen-chuan earthquake of 12 May 2008, Sichuan, People’s Republic of China. GSA Today, 2008, 18(7): 4?11

[25]Cundall P A, Strack O D L. A discrete numerical mode for granular assemblies. Géotechnique, 1979, 29(1): 47?65

[26]Dean S L, Morgan J K, Fournier T. Geometries of frontal fold and thrust belts: insights from discrete element simulations. Journal of Structural Geology, 2013, 53(8): 43?53

[27]Morgan J K, McGovern P J. Discrete element simu-lations of gravitational volcanic deformation: 1. de-formation structures and geometries. Journal of Geo-physical Research: Solid Earth, 2005, 110(B5): 1671? 1677

[28]Naylor M, Sinclair H D, Willett S, et al. A discrete element model for orogenesis and accretionary wed-ge growth. Journal of Geophysical Research: Solid Earth, 2005, 110(B12): 1?16

[29]周易, 于福生, 劉志娜. 分層伸展疊加變形二維離散元模擬. 大地構(gòu)造與成礦學(xué), 2019, 43(2): 213?225

[30]王一丹, 于福生, 劉志娜, 等. 板塊俯沖變形過(guò)程二維離散元模擬——對(duì)東海陸架盆地成因啟示. 海洋地質(zhì)與第四紀(jì)地質(zhì), 2019, 39(5): 166?176

[31]李長(zhǎng)圣. 基于離散元的褶皺沖斷帶構(gòu)造變形定量分析與模擬[D]. 南京: 南京大學(xué), 2019

[32]He R, Shang X, Yu C, et al. A unified map of Moho depth andpsratio of continental China by receiver function analysis. Geophysical Journal International, 2014, 199(3): 1910?1918

[33]Wang W, Wu J, Fang L, et al. Crustal thickness and Poisson’s ratio in southwest China based on data from dense seismic arrays. Journal of Geophysical Research: Solid Earth, 2017, 122(9): 7219?7235

[34]高天揚(yáng), 丁志峰, 王興臣, 等. 利用接收函數(shù)、面波頻散和ZH振幅比聯(lián)合反演青藏高原東南緣地殼結(jié)構(gòu)及其動(dòng)力學(xué)意義. 地球物理學(xué)報(bào), 2021, 64(6): 1885?1906

[35]朱濤, 詹艷, 孫翔宇, 等. 四川龍門(mén)山斷裂帶高精度地殼/巖石圈黏度結(jié)構(gòu)及其動(dòng)力學(xué)意義. 地球物理學(xué)報(bào), 2020, 63(1): 196?209

[36]孫玉軍, 董樹(shù)文, 范桃園, 等. 中國(guó)大陸及鄰區(qū)巖石圈三維流變結(jié)構(gòu). 地球物理學(xué)報(bào), 2013, 56(9): 2936?2946

[37]Diao F, Wang R, Wang Y, et al. Fault behaviorand lower crustal rheology inferred from the first seven years of postseismic GPS data after the 2008 Wen-chuan earthquake. Earth and Planetary Science Letters, 2018, 495: 202?212

[38]Egholm D L, Sandiford M, Clausen O R, et al. A new strategy for discrete element numerical models: 2. sandbox applications. Journal of Geophysical Re-search: Solid Earth, 2007, 112: B05204

[39]Liu Z, Koyi H A. The impact of a weak horizon on kinematics and internal deformation of a failure mass using discrete element method. Tectonophysics, 2013, 586: 95?111

[40]李維波, 李江海, 王洪浩, 等. 庫(kù)車(chē)前陸沖斷帶克拉蘇構(gòu)造帶變形影響因素分析——基于離散元數(shù)值模擬研究. 大地構(gòu)造與成礦學(xué), 2017, 41(6): 1001? 1010

[41]劉暢, 陳曉雪, 張文, 等. PFC數(shù)值模擬中平行粘結(jié)細(xì)觀參數(shù)標(biāo)定過(guò)程研究. 價(jià)值工程, 2017, 36(26): 204?207

[42]Che?ry J, Zoback M D, Hasani R. An integrated me-chanical model of the San Andreas fault in central and northern California. Journal of Geophysical Research: Solid Earth, 2001, 106(B10): 22051?22066

[43]He J K, Lu S J, Wang W M. Thre-dimensional mecha-nical modeling of the GPS velocity field around the northeastern Tibetan plateau and surounding regions. Tectonophysics, 2013, 584: 257?266

[44]Zhu S B, Zhang P Z. FEM simulation of interseismic and coseismic deformation asociated with the 2008 Wenchuan Earthquake. Tectonophysics, 2013, 584: 64?80

[45]尹力, 羅綱. 有限元數(shù)值模擬龍門(mén)山斷裂帶地震循環(huán)的地殼變形演化. 地球物理學(xué)報(bào), 2018, 61(4): 1238?1257

[46]Finch E, Hardy S, Gawthorpe R. Discrete element modelling of contractional fault-propagation folding above rigid basement fault blocks. Basin Research, 2004, 25(4): 515?528

[47]Ai J, Langston P A, Yu H S. Discrete element model-ling of material non-coaxiality in simple shear flows. International Journal for Numerical and Analytical Methods in Geomechanics, 2014, 38: 615?635

[48]Zhang Y, Teng J, Wang Q, et al. Density structure and isostatic state of the crust in the Longmenshan and adjacent areas. Tectonophysics, 2014, 619/620(6): 51?57

[49]Zhu S, Zhang P. FEM simulation of interseismic and coseismic deformation associated with the 2008 Wen-chuan Earthquake. Tectonophysics, 2013, 584: 64?80

[50]Liang S, Gan W, Shen C, et al. Three-dimensional velocity field of present-day crustal motion of the Tibetan Plateau derived from GPS measurements. Journal of Geophysical Research: Solid Earth, 2013, 118(10): 5722?5732

[51]張培震, 鄧起東, 張國(guó)民, 等. 中國(guó)大陸的強(qiáng)震活動(dòng)與活動(dòng)地塊. 中國(guó)科學(xué): D輯, 2003, 33(4): 12?20

[52]Ramberg H. Gravity, deformation and the Earth’s crust in theory, experiments and geological applica-tion. 2nd edition. Kolkata: Academic Press, 1981

[53]Diao F, Wang R, Wang Y, et al. Fault behavior and lower crustal rheology inferred from the first seven years of postseismic GPS data after the 2008 Wen-chuan earthquake. Earth and Planetary Science Let-ters, 2018, 495: 202?212

[54]Huang M H, Bürgmann R, Freed A M. Probing the lithospheric rheology across the eastern margin of the Tibetan Plateau. Earth and Planetary Science Letters, 2014, 396: 88?96

[55]Hao M, Li Y, Zhuang W. Crustal movement and strain distribution in East Asia revealed by GPS observations. Scientific Reports, 2019, 9(1): 16797

[56]瞿偉, 高源, 陳海祿, 等. 利用 GPS 高精度監(jiān)測(cè)數(shù)據(jù)開(kāi)展青藏高原現(xiàn)今地殼運(yùn)動(dòng)與形變特征研究進(jìn)展. 地球科學(xué)與環(huán)境學(xué)報(bào), 2021, 43(1): 182?204

[57]李承濤, 蘇小寧, 孟國(guó)杰. 巴顏喀拉塊體東北緣GPS應(yīng)變率空間分布特征及其與 2017 年九寨溝Ms 7.0地震的關(guān)系. 地震, 2018, 38(2): 37?50

[58]徐錫偉, 聞學(xué)澤, 陳桂華, 等. 巴顏喀拉地塊東部龍日壩斷裂帶的發(fā)現(xiàn)及其大地構(gòu)造意義. 中國(guó)科學(xué):地球科學(xué), 2008, 38(5): 529?543

[59]陳長(zhǎng)云, 任金衛(wèi), 孟國(guó)杰, 等. 巴顏喀拉塊體東部活動(dòng)塊體的劃分、形變特征及構(gòu)造意義. 地球物理學(xué)報(bào), 2013, 56(12): 4125?4141

[60]Wang Z, Wang X, Huang R, et al. Structural hetero-geneities in Southeast Tibet: implications for regio-nal flow in the lower crust and upper mantle. Inter-national Journal of Geophysics, 2012, 315: 1?12

Main Controlling Factors of Uplift Deformation of Longmenshan Structural Belt: Insight from Discrete Element Method

WANG Ying1,2, LI Jianghai1,2,?, MA Changming1,2, SONG Juechen1,2

1. Key Laboratory Orogenic Belts and Crustal Evolution (MOE), School of Earth and Space Sciences, Peking University, Beijing 100871; 2. Institute of Oil and Gas, Peking University, Beijing 100871; ? Corresponding author, E-mail: jhli@pku.edu.cn

In order to explore the main controlling factors of uplift deformation of Longmenshan structural belt, based on the differences in the properties of lower crust material layer between the east and west sides of Longmenshan structural belt, three groups of PFC2D discrete element numerical simulation are carried out to realize quantitative analysis.The experimental deformation results and the model particle motion vector map show that under the condition of no obvious difference in the material properties of the lower crust, the existence of plate collision and compression stress and crustal thickness difference will not form a huge topographic elevation difference in the Longmenshan structural belt. When there are obvious differences in the viscosity coefficient of the lower crust, the relative value of the particle movement rate of the weak lower crust material layer is 1.5?2.94 m/s, and the average movement rate is 1.62 m/s, which is about 54 times of the average movement rate of the particles of the hard lower crust layer. Uplift deformation occurs in the middle of the model (Longmenshan structural belt), with a vertical influence range of 94.74% and a uplift amplitude of 19.85%. The particles of the middle crust and upper crust overlying the weak lower crust have a large upward velocity component, and the upward trend of the material layer of the upper crust is obvious. There is a 20km thickness difference between Bayankala block and the crust of Sichuan Basin, which increases the uplift amplitude of Longmenshan structural belt from 14.79% to 19.85%. Based on the comprehensive analysis of three discrete element simulation experiments, it is concluded that the viscosity difference between the material layer of the lower crust of Bayan Kara block and the material layer of the underground block of Sichuan Basin is the most key control factor for the vertical uplift deformation of Longmenshan structural belt. On the premise that there are obvious differences in the viscosity structure of the lower crust, the crustal thickness differences between the Bayan Kara block and the Sichuan Basin significantly promote the vertical thrust uplift amplitude of the Longmenshan structural belt.

Longmenshan structural belt; discrete element numerical simulation; heterogeneity of the lower crust; crustal thickness; main controlling factors

10.13209/j.0479-8023.2022.056

2021–07–15;

2021–10–15

國(guó)家科技重大專(zhuān)項(xiàng)(2016ZX05033002, 2016ZX05033001)資助