地球物理学报  2019, Vol. 62 Issue (7): 2466-2476   PDF    
褶皱冲断带时空演化动力学过程的数值模拟
周龙寿1,2,3, 刘旭耀4, 胡才博1,2, 孟秋1,2, 张怀1,2, 石耀霖1,2     
1. 中国科学院大学地球与行星科学学院, 北京 100049;
2. 中国科学院计算地球动力学重点实验室, 北京 100049;
3. 中国地震局地壳应力研究所, 北京 100085;
4. 中国地质大学(武汉)工程学院, 武汉 430074
摘要:利用库仑临界锥角理论和沙箱物理模拟进行褶皱冲断带的研究时,通常忽略了楔形体介质的内聚力.岩石力学实验结果表明,岩石内聚力通常在几到几十兆帕范围内变化.楔形体介质强度的变化是否会影响褶皱冲断带的时空演化?针对这一问题,我们建立了岩石内聚力分别是10 MPa、20 MPa和50 MPa的3个二维弹塑性有限元模型,模型包含了楔形体介质的弹塑性材料非线性和底部滑脱带的接触非线性.该模型考虑了不同介质强度、底部滑脱带摩擦、重力和边界构造加载的影响,更为符合实际.计算结果表明,岩石内聚力为10 MPa时,楔形体内的断坡首先在楔形体近端产生,从近端依次向远端发展;岩石内聚力是20 MPa时,断坡开始在楔形体近端产生,随后在远端也开始形成,最后由两端向中间汇聚;岩石内聚力是50 MPa时,断坡先从楔形体远端形成,从远端向近端依次发展.我们还讨论了底部滑脱层倾角对褶皱冲断带演化的影响,结果表明较低的底部滑脱面倾角易产生由近端向远端演化的样式,中等滑脱面倾角易产生两端向中间演化的样式,较高滑脱面倾角易产生由远端向近端演化的样式.我们得到了三种不同的褶皱冲断带时空演化的模式,其结果可以用来解释青藏高原东北缘依次向北东方向发展的海原-六盘山断层、天景山断层、烟筒山断层系统的时空演化.
关键词: 楔形体介质      底部滑脱层      弹塑性模型      褶皱冲断带      青藏高原东北缘     
Numerical simulation of the dynamic process of spatial and temporal evolution for fold-thrust belts
ZHOU LongShou1,2,3, LIU XuYao4, HU CaiBo1,2, MENG Qiu1,2, ZHANG Huai1,2, SHI YaoLin1,2     
1. School of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing 100049, China;
2. Key Laboratory of Computational Geodynamics, Chinese Academy of Sciences, Beijing 100049, China;
3. Institute of Crustal Dynamics, China Earthquake Administration, Beijing 100085, China;
4. Engineering College of China University of Geosciences(Wuhan), Wuhan 430074, China
Abstract: The effect of cohesion is usually neglected in the study of deformation of fold and thrust belts by using Coulomb critical wedge theory and sandbox physical simulation. Rock mechanic experiments show that the cohesion of rocks varies from several to several dozen MPa. Will the change of the strength of a wedge medium affect the temporal and spatial evolution of fold and thrust belts? To address this issue, we construct three two-dimensional elastoplastic finite element models containing rock cohesion of 10 MPa, 20 MPa and 50 MPa, respectively. The models include nonlinear elastoplastic material of a wedge-shaped medium and contact nonlinear nature of the bottom detachment zone. The effects of different medium strength, bottom detachment zone, gravity, and boundary tectonic loading are considered in the models, which can mimic realistic geology. The results show that when the rock cohesion is 10 MPa, a thrust ramp in the wedge first appears at the near end of the wedge and develops from the near end to the distant end in turn. When the rock cohesion is 20 MPa, the thrust ramp begins to appear at the near end of the wedge body, then begins to form at the distant end, and finally converges from the two ends to the middle. When the rock cohesion is 50 MPa, the ramp begins to form from the distant end of the wedge, and develops from the distant end to the near end. We also examine the influence of the dip angle of the basal detachment on the evolution of fold-thrust belts. The results show that the lower dip angle of the basal detachment tends to evolve from the near end to the distant, the middle dip angle of the basal detachment tends to evolve from the two ends to the middle, and the higher dip angle of the basal detachment tends to evolve from the distant to the near end. We have obtained three different models of temporal and spatial evolution of fold-thrust belts and examined the influence of wedge-shaped medium strength on the evolution of thrust ramps. The results can be used to explain the temporal and spatial evolution of the Haiyuan-Liupanshan, Tianjingshan and Yantongshan fault systems in the northeastern margin of the Tibet plateau, which are developing northward and eastward in turn.
Keywords: Wedge medium    Basal detachment    Elastoplastic model    Fold-and-thrust belt    Northeastern margin of Tibetan Plateau    
0 引言

褶皱冲断带是一种在挤压构造加载和底部滑脱层活化的情况下,楔形体内部产生的以褶皱和逆冲断层为主的强烈变形地质构造带,主要发育在汇聚型板块边界,如陆-陆碰撞或陆内缩短、前陆造山带等.褶皱冲断带中形成的一系列逆冲断层称之为断坡,断坡在时间和空间分布上都有一定的规律性.

褶皱冲断带这一概念首先应用于加拿大的落基山(Price and Mountjoy, 1971),20世纪70、80年代随着油气勘探的重大突破,褶皱冲断带引起越来越多的关注及深入研究(Chapple,1978Davis et al., 1983; Davis and Engelder, 1985).根据这些勘探成果,可以概括出褶皱冲断带的3个基本特征:1)褶皱冲断带内后陆位置强烈增厚,可以吸收上百公里的缩短,到前缘却只有几公里的变形;2)褶皱冲断带存在一个向后陆倾斜的底部滑脱面或者叫拆离断层面,该面以下的介质变形微弱,该面以上的介质褶皱和逆冲推覆构造强烈;3)褶皱冲断带垂直剖面中滑脱面的上部地层一般为楔形形态,底部滑脱面向后陆倾斜,顶面坡角主要向前陆倾斜,剖面整体形态向前陆收敛.

关于褶皱冲断带形成的力学机制研究有3种主流方法:1)库仑临界锥角理论;2)沙箱物理模拟;3)数值模拟.库仑临界锥角理论可以解释一系列向前发展的褶皱冲断带和增生楔的发展演化,但通常忽略楔形体内聚力并假设楔形体内部和底部滑脱带介质均匀(Chapple, 1978; Davis et al., 1983; Dahlen et al., 1984; Dahlen, 1990; Bilotti and Shaw, 2005; Cubas et al., 2008);沙箱物理模拟沉积盖层普遍采用粒径为200~400 μm的松散石英砂,内聚力接近于0(Storti and McClay, 1995; Storti et al., 2000; Costa and Vendeville, 2002; Agarwal and Agrawal, 2002; Lujan et al., 2003; Bonnet et al., 2007; Bose et al., 2009; Zhou et al., 2007; Malavieille, 2010; Mary et al., 2013);利用数值模拟方法研究褶皱冲断带是一种非常有效的手段,可以定量、形象地揭示褶皱冲断带个性化的几何学、运动学、动力学等各方面特征.前人的数值模拟模型主要分为离散元(Hardy and Finch, 2005)、有限差分(Ruh et al., 2013)和有限元(Beaumont et al., 2001; Panian and Wiltschko, 2004, 2007; Ruh et al., 2012; Fillon et al., 2013)3种.

前人已经探讨了影响褶皱冲断带时空演化的各种因素.

首先是底部滑脱层深度的影响.若滑脱层存在于地壳的脆韧转换带以下,则由于温度或者岩石类型不同,底部滑脱层呈现出与压力无关的塑性现象,基于压力的临界库仑楔理论已经不适用;若滑脱层存在脆韧转换带以上的上地壳中,则造山过程适用临界库仑楔理论(Davis et al., 1983; Dahlen, 1990; Liu et al., 1992).第二,孔隙流体压力在断坡形成机制中起着重要作用,需要注意陆上和海面下的楔形体形成过程中孔隙流体压力的不同影响(Hubber and Rubey, 1959).

底部滑脱带的剪切强度和倾角对褶皱冲断带形态的影响很大,剪切强度的高低直接影响褶皱冲断带演化过程是否符合库仑楔理论,滑脱层倾角大小则影响了楔形体内断坡的汇聚形态(Davis et al., 1983; Dahlen, 1990; Liu et al., 1992; Smit et al., 2003).底部滑脱层剪切强度较低,褶皱冲断带由两端向中间汇聚,底部滑脱层比较坚硬时,向前发展冲断带占优势,褶皱发展不对称,呈旋转状分布(Costa and Vendeville, 2002),单一摩擦基底的褶皱冲断带形成逆冲系统,而黏性多滑脱层的结构演变取决于基底和顺序滑脱层之间的强度关系,单一摩擦基底和多滑脱层的结果和莫尔—库仑临界楔理论一致(Ruh et al., 2012).Zhang等(2018)利用极限分析理论,分别采用解析解和数值模拟的方法,考虑了深部和浅部两个滑脱层的存在对青藏高原东缘龙门山地区和四川盆地内一系列逆冲断层系统时空演化的影响.

短时间内同步发生的沉积和剥蚀能够保持临界楔的形态,前陆盆地的演化取决于被侵蚀和沉积的楔形体数量,前陆盆地的演化和内部结构由楔形体形成机制所控制(Storti and McClay, 1995Mugnier et al., 1997Bonnet et al., 2007).底部沉积物的陆向快速加载可能会影响底部滑脱层的流变特性,与增生楔底部接触的黏弹性地层也有可能会在板间大地震中破裂(Gutscher et al., 1998, 2001).地形起伏对褶皱冲断带发展也有很大的影响,在没有侵蚀和沉积情况下,地表坡度趋于均匀,与临界楔理论预测相同.侵蚀或沉积则使地形剖面不规则,冲断带无序发展(Marques and Cobbold, 2002).

数值模拟方法的优点是定解问题的方程清楚,研究区域可以很大,边界条件根据实际地质边界加载给定,本构方程较为多样,材料参数由实验确定,操作简单,可重复.褶皱冲断带的数值模拟一般考虑了1)重力;2)构造加载;3)地形;4)滑脱带倾角;5)介质不均匀;6)内摩擦角;7)滑脱带的强度等因素的影响,但介质不均匀性不是必须的(Panian and Wiltschko, 2007).本文在前人工作基础上,考虑了楔形体重力和构造加载,建立楔形体内聚力分别为10 MPa、20 MPa和50 MPa的带有底部滑脱接触带的二维弹塑性有限元模型,研究了上部楔形体介质强度不同时增生楔内断坡的形成和演化过程,这个研究结果有助于解释青藏高原东北缘依次向北东方向发展的海原—六盘山断层、天景山断层、烟筒山断层的时空演化.

1 有限元模型

本文的有限元模型包含两个非线性:材料非线性和接触非线性.与一般的线弹性有限元模型相比,本文采用的模型考虑了上部楔形体介质强度以及底部滑脱带摩擦的影响,更符合实际,但求解难度比较大,可能存在计算不稳定的情况.

1.1 理论介绍

模型中采用理想弹塑性本构关系,使用Drucker-Prager破裂准则(Panian and Wiltschko, 2004, 2007)如下所示:

(1)

式中,I1是应力张量σij的第一不变量,I1=σkkJ2是偏应力张量的第二不变量,J2=α, k是材料参数.

Drucker-Prager屈服面在主应力空间中是一个圆锥面,它在σ1+σ2+σ3=0的平面(称为π平面)上的截线是一个圆.为确定参数αk与岩土工程计算中常用的内聚力c和摩擦角φ之间的关系,需要将Drucker-Prager屈服条件式(1)与如下的库仑屈服条件相比较,即

(2)

如果不规定主应力的大小次序σ1σ2σ3,采用对称开拓法,库仑屈服条件在主应力空间中是一个六棱锥.它在π平面内的截线是一个不规则的六边形.使Drucker-Prager屈服面的锥点和库仑梭锥的顶点重合,Drucker-Prager屈服面在π平面内的截线称为Drucker-Prager圆,当它与库仑六边形的外顶点重合时,可得

(3)

当它与库仑六边形的内顶点重合时,可得

(4)

Drucker-Prager屈服条件考虑了围压(静岩压力)对屈服特性的影响,并且能反映剪切引起膨胀(扩容)的性质,在模拟岩石材料的弹塑性性质时,这种屈服条件得到广泛的应用.

等效应变增量dε(王平和崔建忠,2006)由下式计算得出:

(5)

等效塑性应变增量dεp的计算公式为

(6)

式中下标p表示塑性变形引起的塑性应变.总的等效塑性应变εp通过不同载荷步的积分运算得到:

(7)

1.2 几何模型

在Comsol Multiphysics软件中建立二维弹塑性有限元模型,垂直剖面如图 1所示.

图 1 楔形体几何模型及边界条件示意图 Fig. 1 Geometric model and boundary conditions of wedge

几何模型分为上下两个部分.上部为楔形体,下部为几乎不变形的下覆地层.上部楔形体水平方向长100 km,左端为厚端,垂直厚度10 km.楔形体地表倾角2°,底部滑脱带倾角2°,两者倾向相反.下覆地层水平长度为130 km,左端为薄端,垂直厚度10 km.楔形体和下覆地层之间通过一个底部滑脱层连接.如图 1所示,上部楔形体的左端为边界构造加载端,称之为近端,几何上称为厚端;楔形体的右端的水平位移约束为零,称之为远端,几何上是薄端.

1.3 材料参数

上下两部分的介质均采用弹塑性材料,应用Drucker-Prager破裂准则及非关联流动法则(Panian and Wiltschko, 2004, 2007).楔形体与下覆地层之间采用了没有厚度的库仑-莫尔摩擦接触面,摩擦系数为0.175(Panian and Wiltschko, 2007).这种接触面是高度非线性的,因为摩擦接触力和接触状态事先是未知的,需要反复迭代才能确定.通常有两种方法来确定摩擦接触力和摩擦滑动的位移,一种是使用罚函数法,需要引入刚度极大的弹簧,优点是不引入新的求解变量,缺点是弹簧刚度的选取具有人为性.另一种是采用拉格朗日不连续变形分析方法(LDDA),在数学上引入拉格朗日乘子这个新的变量,需要补充位移的相关约束条件才能求解(Cai et al., 2000蔡永恩等,2004).3个模型的具体材料参数见表 1.

表 1 模型材料参数 Table 1 Model material parameters
1.4 边界条件

模拟分为两步.第一步只考虑岩石自重,不考虑边界加载.如图 1所示.模型上下两个部分的左右边界,都是位移滚筒约束,法向位移约束为零,下覆地层的底部完全固定.上部楔形体表面自由.重力从100 N·m-3逐步施加,一个加载步施加200 N·m-3,共加载123步,最终达到24500 N·m-3.第一步最终的位移场和应力场作为第二步构造加载的初始场.模拟的第二步,在第一步重力全部施加后,从上部楔形体的左侧边界进行水平构造加载.一个加载步加载10 m向右的水平位移增量,一共加载600步,水平构造加载共计6000 m.其他边界条件与第一步相同.

重力是体力,在自然界中直接一次性作用于地球介质的每一点.对于数值模型而言,如果考虑重力的影响,重力的施加有两种方式:一种是一次性作为体力施加到研究区域,但必须是线性模型,比如弹性模型,或者黏弹性模型,即一次性加载与分步加载对于线性模型而言,结果是一样的;另一种是分步施加重力,对应于非线性问题,比如弹塑性模型、接触非线性模型.本文采用的是弹塑性模型和接触非线性模型,必须采用分步加载重力的方式施加重力,否则容易产生较大的塑性变形或者接触力计算不准,从而导致数值计算不稳定或者失效.

2 模拟结果 2.1 模型1(内聚力为10 MPa)

图 2是5个不同构造加载阶段上部楔形体内等效塑性应变增量的结果.等效塑性应变是衡量介质达到塑性屈服之后塑性变形大小的重要指标.当左侧构造加载增加到300 m时,楔形体内开始产生较明显的塑性变形.以300 m为第一个断坡演化的起始阶段,将整个断坡时空演化分为5个构造阶段:300~1300 m,1300~2300 m,2300~3300 m,3300~4300 m,4300~5300 m.最后看到,当构造加载到5300 m时,模型1的上部楔形体内共形成了5个相对完整的从厚端(近端)向薄端(远端)依次演化的断坡,相邻断坡之间还形成较为明显的塑性变形集中带(褶皱),它们共同构成了楔形体内的褶皱冲断带.

图 2 楔形体中不同加载阶段内的等效塑性应变增量分布图 (a) 300~1300 m;(b) 1300~2300 m;(c) 2300~3300 m;(d) 3300~4300 m;(e) 4300~5300 m. Fig. 2 Distribution of equivalent plastic strain increment in wedge under different loading stages

图 2中等效塑性应变增量的分布可以很清楚地看到楔形体内褶皱冲断带和断坡的时空演化过程.图 2中的色标都是0~0.1(下同,等效塑性应变没有量纲),等效塑性应变会在某些特定位置显著增加并成条带状分布,这些条带状分布的区域就是潜在的褶皱冲断带断坡和褶皱的位置.断坡和褶皱的具体位置分布如下:

图 2a清楚地显示第一条断坡的形成位置,对应于加载的第一阶段(300~1300 m).第一条断坡从楔形体最左端底部滑脱层开始激发,逐渐发展到楔形体上表面,断坡倾角约30°.

图 2b显示了第二条断坡和相邻的褶皱的形成位置,对应于构造加载的第二阶段(1300~2300 m).与第一条断坡(图 2a)相比,第二条断坡也是从底部滑脱带开始激发,以几乎相同的倾向发展到楔形体顶部.在构造加载的第二阶段(1300~2300 m),第一条断坡和第二条断坡之间,存在着一条较为明显的等效塑性应变集中带,连接第一条断坡的顶部和第二条断坡的底部,对应于地质上的褶皱变形.可以看出,第二条断坡内的等效塑性应变增量比左侧相邻的褶皱内的要大一些,呈现优先发展的趋势.前两个阶段形成的断坡近似平行,两者在地表的水平距离约为30 km,约占整个上部楔形体水平长度的3/10.

图 2c显示了第三条断坡和相邻的褶皱的形成位置,对应于构造加载的第三阶段(2300~3300 m).与第二阶段类似,第三条断坡内的等效塑性应变比左侧相邻的褶皱内的要大一些,并且断坡也是从底部滑脱带开始激发,以跟前两条断坡几乎相同的倾向发展到楔形体顶部.第三阶段内形成的褶皱与第二阶段内形成的断坡也在楔形体顶部汇聚.三个阶段形成的三个断坡近乎平行,依次从厚端向薄端传播.

图 2(de)显示了第4个加载阶段和第5个加载阶段内形成的褶皱和断坡.结果与第二个和第三个加载阶段类似.

总体来看,从加载300 m开始至加载到5300 m,在水平长度为100 km的楔形体内形成了5个断坡.从300 m到5300 m大致可分为5个加载阶段,第2~5个加载阶段,不仅产生了断坡,还相应产生了褶皱,共同形成了褶皱冲断带.褶皱和断坡都是从底部滑脱带产生,逐渐以相反倾向延伸到楔形体顶部.后一阶段的褶皱和前一阶段的断坡在楔形体顶部汇聚.5条断坡几乎平行分布,从厚端向薄端依次演化,相邻两条断坡的水平距离随加载而逐渐减小,依次为:30 km,20 km,15 km,10 km.

图 3给出了5个不同加载阶段垂直位移增量的分布.从图 3a图 3e可以看出,在5个不同的构造加载阶段,在褶皱和断坡相夹的部分,存在明显的垂直向上的位移分量,表明该部分介质在增厚.除了褶皱和断坡相夹的部分存在明显垂直位移之外,楔形体其余部分垂直位移并不明显,也就是说楔形体内增厚的部分集中在褶皱和断坡相夹的部分.对每一个褶皱和断坡相夹的部分,在其中间部位增厚最多,由中间向两边增厚逐渐减小.垂直增厚从楔形体的厚端向薄端依次减小,地表增厚从第一个加载阶段的800 m降低到第5个加载阶段的100 m.垂直增厚的部位并不是出现在断坡出露地表处,而是出现在相邻两个断坡出露地表处的中间,这是由于本文采用的连续变形的上部楔形体模型,断坡的塑性变形虽然非常显著,但仍然是连续变形,并没有发生错动,还有一个原因本文没有考虑地表的剥蚀和沉积作用,不考虑地形起伏引起的地表物质运移.

图 3 楔形体中不同加载阶段内垂直位移增量分布图 (a) 300~1300 m;(b) 1300~2300 m;(c) 2300~3300 m;(d) 3300~4300 m;(e) 4300~5300 m. Fig. 3 Distribution of vertical displacement increment in wedge under different loading stages

图 4给出了5个不同加载阶段水平位移增量的分布.从图 4可以看出,水平位移从厚端开始出现,随着边界加载位移的不断增加逐渐向薄端发展.和垂直位移只出现在褶皱和断坡相夹部分不同,水平位移在整个楔形体内均有发生.楔形体内水平位移基本按照断坡和褶皱的发生位置分块分布,随着边界加载位移的加大,各块体内的水平位移也不断加大.总体而言,离加载端越近,水平位移越大;离加载端越远,水平位移越小.

图 4 楔形体中不同加载阶段内水平位移增量分布图 (a) 300~1300 m;(b) 1300~2300 m;(c) 2300~3300 m;(d) 3300~4300 m;(5) 4300~5300 m. Fig. 4 Distribution of horizontal displacement increment in wedge under different loading stages
2.2 模型2(内聚力为20 MPa)

模型2与模型1相比,仅仅是上部楔形体的内聚力从10 MPa增加到20 MPa,其余材料参数不变,几何模型和边界条件均不发生变化.图 5显示了模型2内聚力为20 MPa时楔形体内不同加载阶段的等效塑性应变增量.

图 5 模型2楔形体中不同加载阶段等效应变增量分布图 (a) 300~1300 m; (b) 1300~2300 m; (c) 2300~3300 m; (d) 3300~4300 m; (e) 4300~5300 m. Fig. 5 Distribution of equivalent strain increment at different loading stages in model 2 wedge

图 5a可以看出,加载起始阶段(300~1300 m),楔形体内断坡出现位置、形态和模型1楔形体内聚力为10 MPa时类似,也是由底部滑脱层开始激发,向楔形体上表面开始发展,到距厚端20 km处,断坡贯通到楔形体上表面.

图 5b图 4不同,加载至1300~2300 m时,楔形体内除了厚端(近端)继续出现与前一个断坡汇聚的褶皱和断坡外,最远端(薄端)也出现了断坡和褶皱,断坡由楔形体的上顶面和坡底面同时激发,最后上下贯通,形成一个完整的断坡,之后这个断坡向近端发展出倾向相反的褶皱,褶皱与断坡共轭出现,与断坡汇聚于坡底.

图 5c显示了构造加载为2300~3300 m时的情况.随着构造加载距离加大,厚端形成的第一条断坡等效塑性应变增量减小,但远端断坡和褶皱的塑性应变增量增大.同时楔形体中间又发展出3条断坡,褶皱也随之共轭出现,分别与前一条断坡汇聚于坡面,但靠近远端的第3条断坡不是很明显,塑性应变增量较小.靠近远端断坡应变较大,图 5b中远端出现的第一条断坡塑性应变增大.

图 5d图 5c类似,只是靠近远端的第二条断坡在坡顶处的塑性应变增量增大.图 5e显示随着构造加载距离至4300~5300 m,楔形体内最终形成5个褶皱及其相邻的褶皱.

总体来说,内聚力为20 MPa时,楔形体内的褶皱冲断带不是从近端向远端一个方向发展,而是由两端向中间汇聚,相邻的褶皱和断坡在楔形体坡顶和坡底相互连通,模型2中断坡演化的空间顺序与模型1(内聚力为10 MPa)从近端到远端单方向演化的空间顺序不同.

2.3 模型3(内聚力为50 MPa)

图 6显示的是模型3楔形体内聚力为50 MPa时楔形体内部不同加载阶段等效塑性应变增量.结果显示,褶皱和断坡从楔形体远端(薄端)开始发育.加载至300~1300 m时(图 8a),前端发育出第一条褶皱冲断带,由楔形体上表面开始激发,延伸至坡底后,又由坡底向坡顶延伸,长度大约10 km.随着近端(厚端)的不断加载,褶皱冲断带不断由薄端向厚端发育,等效塑性应变增量的集中条带比较明显.到加载至5300 m时,楔形体内部由远端向近端依次发育出了3条完整的褶皱冲断带,同时近端出现了1条塑性应变增量较小的断坡.模型3(楔形体内聚力为50 MPa)和模型1(内聚力为10 MPa)的区别在于褶皱冲断带演化方向正好相反,前者楔形体内褶皱冲断带基本是从远端向近端发育,后者褶皱冲断带基本是从近端向远端发育,而模型2(楔形体内聚力为20 MPa)楔形体内褶皱冲断带基本是从两端向中间发育.

图 6 模型3楔形体内聚力为50 MPa时楔形体内不同加载阶段等效塑性应变增量 (a) 30 m~1300 m; (b)1300~2300 m; (c) 2300~3300 m; (d)3300~4300 m; (e)4300~5300 m. Fig. 6 Equivalent plastic strain increment at different loading stages in model 3 wedge with cohesion of 50 MPa
图 8 上部楔形体内聚力为20 MPa,不同滑脱面倾角β对褶皱冲断带演化的影响(其余参数与模型2相同) (a) β=1°,左侧边界加载4 km;(b) β=2°,左侧边界加载4 km;(c) β=4°,左侧边界加载4 km;(d) β=8°,左侧边界加载1.05 km. Fig. 8 The influence of different dip angles β of decollements on the evolution of fold-thrust belts when the cohesion in the upper wedge is 20 MPa (the other parameters are the same as those in Model 2) (a) β=1°, loading 4 km on the left boundary; (b) β=2°, loading 4 km on the left boundary; (c) β=4°, loading 4 km on the left boundary; (d) β=8°, loading 1.05 km on the left boundary.
3 讨论 3.1 褶皱冲断带时空演化的力学原理

楔形体在构造应力作用下的力学原理图如图 7所示,根据水平方向力的平衡条件,P=R+fPR,左侧的水平推力要大于右侧的水平力.但是从应力的定义来看,对于平面应变问题(模型为单位厚度),σleft=P/H0σright=R/H1H0H1.当楔形体的强度较小,内聚力较小,内摩擦角较小,滑脱层的摩擦系数较小时,应力很难从近端直接传递到远端,越靠近近端,H1越接近H0,此时σleftσright,这种情况最容易从近端产生塑性变形集中带,并逐渐向远端扩展生长,对应本文模型一的情形.

图 7 褶皱冲断带变形的力学机制原理图 P为近端水平推力,对应底层厚度为H0R为底层内水平力,对应地层厚度为H1f为底部滑脱面的摩擦力的水平分量;α、β分别为楔形体的顶面倾角和底面倾角. Fig. 7 Principle diagram of mechanics mechanism of deformation in fold-thrust zone P is the proximal horizontal thrust, the corresponding stratum thickness is H0, R is the horizontal force in the stratum, the corresponding stratum thickness is H1, f is the horizontal component of the friction force on the bottom decollement; α and β are the top and bottom dip angles of the wedge respectively.

当楔形体的强度较大,内聚力较大,内摩擦角较小,滑脱层的摩擦系数较小时,应力比较容易从近端直接传递到远端,远端的水平力小于近端的水平推力P,但是远端的H1H0小,此时σleftσright,这种情况最容易从远端产生塑性变形集中带,并逐渐向近端扩展生长,对应本文模型三的情形.

当模型的强度适中时,则可能出现σleft=σright,近端和远端同时出现塑性变形集中带,并向模型中间扩展生长,对应本文模型二的情形.

本文中的几何模型分上下两部分,上部是楔形体,楔形体近端较厚,为10 km,远端较薄,约3 km,代表比较松散的沉积层,在长期构造加载的作用可以发生较大的变形,发育一系列的褶皱冲断带.下部为比较致密的基岩层,变形较小,本文通过较大的内聚力和内摩擦系数以提高其介质强度,在重力及不断增加的边界构造加载的共同作用下,在底部滑脱层的摩擦接触条件下,上部楔形体逐渐产生塑性变形集中条带,这些塑性变形集中条带与地质上的褶皱冲断带较为类似.

Smit等(2003)总结了褶皱冲断带的演化形式存在的4种类型,其中最常见的就是从近端向远端向前发展的褶皱冲断带类型,其特点是依次向前发展的断层为主断层,对应本文模型得到的断坡,同时存在与主断层垂直的向后发展的转换断层,对应于本文模型得到的断坡相邻的褶皱.Panian and Wiltschko(2007)利用弹塑性模型模拟了从近端向远端依次向前发展的褶皱冲断带的时空演化,对应于本文模型1(楔形体内聚力为10 MPa)的结果,本文还给出了当楔形体内聚力较大时(20 MPa和50 MPa)褶皱冲断带的时空演化,结果表明存在三种不同的褶皱冲断带的时空演化模式.

除介质强度外,楔形体顶角倾角也对褶皱冲断带的形成演化有着巨大的影响.为了研究滑脱面倾角β对褶皱冲断带演化的影响,我们在模型2的基础上,在保证上部楔形体内聚力为20 MPa,其余模型模型参数也不变的情况下,仅仅改变滑脱面倾角,建立了如下4个模型,其滑脱面倾角分别为1°, 2°, 4°, 8°.图 8给出了4个模型总的塑性应变分布对比图,其中前三个模型的左侧边界加载为4 km,最后一个模型的左侧边界加载为1.05 km.滑脱面倾角为8°时,因为数值模型弹塑性材料的非线性和滑脱面的接触非线性,模拟计算当左侧边界加载达到1.05 km时计算不收敛.

图 8的结果显示,滑脱面倾角对褶皱冲断带的影响很大.如上文提及的三个模型,β=1°(低倾角)时,褶皱冲断带从近端向远端逐渐发展(图 8a),类似于模型1(低内聚力10 MPa);β=2°(中等倾角)时,就是模型2本身(中等内聚力20 MPa),褶皱冲断带从远近两端向中间逐渐发展(图 8b);β=4°(高倾角)时,褶皱冲断带从远端向近端逐渐发展(图 8c),类似于模型3(高内聚力50 MPa);当滑脱面倾角达到β=8°(高倾角)时,褶皱冲断带也是从远端向近端逐渐发展(图 8c),因为模型的材料非线性和接触非线性,数值计算在左侧边界加载到1.05 km时发散不收敛.

楔形体底部滑脱面的内摩擦角也是褶皱冲断带时空演化的重要影响因素之一.当内摩擦角较小(φb=5°)时,楔形体在构造加载作用下的变形(应变率张量的第二不变量)较大,可以由近端一直传递到远端;随着内摩擦角增大(φb=10°、φb=15°)时,楔形体内褶皱冲断带的变形由近端传向远端的距离也在减小,整个楔形体只有在靠近构造加载的近端至楔形体中部发生变形,中部至远端则没有任何变化(Ruh et al., 2013).

由以上讨论可知,楔形体内褶皱冲断带的时空演化模式由模型的塑性强度、底部滑脱层的摩擦系数、楔形体顶底倾角几何参数等因素共同控制.

3.2 青藏高原东北缘褶皱冲断带形成发展的可能成因

Wang等(2013)通过大量的地质考察、岩石年代学研究了宁夏盆地南部自第三纪中新世-上新世盆地造山及逆冲断层系统时空演化过程,约10 Ma六盘山和南华山因为海原—六盘山逆冲断层的北东向逆冲而快速抬升造山,约5.4 Ma由于天景山断裂带的活化导致天景山的抬升,约2.6 Ma烟筒山断裂带的北东向逆冲导致烟筒山的抬升,这整个过程造成青藏高原东北缘不断向北东生长扩展.本文建立的模型1能够用来解释青藏高原东北缘依次向北东方向发展的海原—六盘山断层、天景山断层、烟筒山断层的时空演化,暗示盆地内10~2.6 Ma的地质时间内浅层沉积层的岩石强度较低,岩石内聚力较小,地壳内的滑脱面倾角较小.

青藏高原东北缘六盘山地区基底下方是沉积物和板块碎片(洋壳和陆壳碎片)所组成的构造增生楔(郭晓玉等,2017);通过对长约90km的深地震反射剖面研究,得出了海原断裂的深部几何形态和两侧精细地壳结构(王海燕等,2012);海原断裂下方波速比有较异常的增大(沈旭章等,2013),说明海原断裂下方地壳中存在松软的物质.除以上文献外,前人关于青藏高原东北缘地壳结构的工作(李文辉等,2017陈九辉等,2005)中,均对海原—六盘山断裂带或包括该断裂在内的更大范围区域研究比较精细,结果表明海原—六盘山断裂区域地壳中存在比较松软的低波速沉积物质.”

我们的数值模拟工作根据实际观测资料进行.由于该区域的天景山断层、烟筒山断层及罗山—牛首山断裂属于同一盆地内的较小断裂,目前对这些断裂的观测精度还比较粗略,上覆地层的最大深度为10 km,地震层析成像结果和高精度地震反射剖面结果对研究区域的盆地沉积层的物性差异没有足够的分辨率,岩石钻探和采样的数据较为缺乏,相关岩石物理实验资料也较少.我们将在该区域地质观测和岩石钻探、采样及力学实验资料更为精细丰富之后,展开强度差异影响的研究.

4 结论

本文通过有限元数值模拟方法研究了褶皱冲断带的形成及演化过程,结果表明上部楔形体的强度和底部滑脱面的倾角是其重要的影响因素,存在三种不同的褶皱冲断带时空样式.较低强度楔形体和较低的底部滑脱面倾角易产生由近端向远端演化的样式,中等强度楔形体和中等滑脱面倾角易产生两端向中间演化的样式,较高强度楔形体和较高滑脱面倾角易产生由远端向近端演化的样式.较低强度楔形体和较低的底部滑脱面倾角的模型可用于解释青藏高原东北缘宁夏盆地南部自第三纪中新世-上新世海原—六盘山、天景山、烟筒山断层系统的时空演化过程.

References
Agarwal K K, Agrawal G K. 2002. Analogue sandbox models of thrust wedges with variable basal frictions. Gondwana Research, 5(3): 641-647. DOI:10.1016/S1342-937X(05)70635-3
Beaumont C, Jamieson R A, Nguyen M H, et al. 2001. Himalayan tectonics explained by extrusion of a low-viscosity crustal channel coupled to focused surface denudation. Nature, 414(6865): 738-742. DOI:10.1038/414738a
Bilotti F, Shaw J H. 2005. Deep-water Niger Delta fold and thrust belt modeled as a critical-taper wedge:The influence of elevated basal fluid pressure on structural styles. AAPG Bulletin, 89(11): 1475-1491. DOI:10.1306/06130505002
Bonnet C, Malavieille J, Mosar J. 2007. Interactions between tectonics, erosion, and sedimentation during the recent evolution of the Alpine orogen:Analogue modeling insights. Tectonics, 26(6): TC6016. DOI:10.1029/2006TC002048
Bose S, Mandal N, Mukhopadhyay D K, et al. 2009. An unstable kinematic state of the Himalayan tectonic wedge:Evidence from experimental thrust-spacing patterns. Journal of Structural Geology, 31(1): 83-91. DOI:10.1016/j.jsg.2008.10.002
Cai Y E, He T, Wang R. 2000. Numerical simulation of dynamic process of the Tangshan earthquake by a new method-LDDA. Pure and Applied Geophysics, 157(11-12): 2083-2104.
Cai Y E, Luo G, Liu J Z, et al. 2004. A numerical method for modeling dynamic process of discontinuous rock foundation-arch dam system. Rock and Soil Mechanics (in Chinese), 25(S2): 361-365.
Chapple W M. 1978. Mechanics of thin-skinned fold-and-thrust belts. Geological Society of America Bulletin, 89(8): 1189-1198. DOI:10.1130/0016-7606(1978)89<1189:MOTFB>2.0.CO;2
Chen J H, Liu Q Y, Ii S C, et al. 2005. Crust and upper mantle S-wave velocity structure across Northeastern Tibetan Plateau and Ordos block. Chinese Journal of Geophysics (in Chinese), 48(2): 333-342.
Costa E, Vendeville B C. 2002. Experimental insights on the geometry and kinematics of fold-and-thrust belts above weak, viscous evaporitic décollement. Journal of Structural Geology, 24(11): 1729-1739. DOI:10.1016/S0191-8141(01)00169-9
Cubas N, Leroy Y M, Maillot B. 2008. Prediction of thrusting sequences in accretionary Wedges. Journal of Geophysical Research:Solid Earth, 113(B12): B12412. DOI:10.1029/2008JB005717
Dahlen F A, Suppe J, Davis D. 1984. Mechanics of fold-and-thrust belts and accretionary wedges cohesive coulomb theory. Journal of Geophysical Research, 89(10): 100087-10101.
Dahlen F A. 1990. Critical taper model of fold-and-thrust belts and accretionary wedges. Annual Review of Earth and Planetary Sciences, 18: 55-99. DOI:10.1146/annurev.ea.18.050190.000415
Davis D, Suppe J, Dahlen F A. 1983. Mechanics of fold-and-thrust belts and accretionary wedges. Journal of Geophysical Research:Solid Earth, 88(B2): 1153-1172. DOI:10.1029/JB088iB02p01153
Davis D M, Engelder T. 1985. The role of salt in fold-and-thrust belts. Tectonophysics, 119(1-4): 67-88. DOI:10.1016/0040-1951(85)90033-2
Fillon C, Huismans R S, Van Der Beek P. 2013. Syntectonic sedimentation effects on the growth of fold-and-thrust belts. Geology, 41(1): 83-86. DOI:10.1130/G33531.1
Guo X Y, Gao R, Gao J R, et al. 2017. Integrated analysis on the deformation of the Liupan Shan fold-thrust belt, NE Tibet, and its tectonic attribution. Chinese Journal of Geophysics (in Chinese), 60(6): 2058-206. DOI:10.6038/cjg20170603
Gutscher M A, Kukowski N, Malavieille J, et al. 1998. Material transfer in accretionary wedges from analysis of a systematic series of analog experiments. Journal of Structural Geology, 20(4): 407-416. DOI:10.1016/S0191-8141(97)00096-5
Gutscher M A, Klaeschen D, Flueh E, et al. 2001. Non-Coulomb wedges, wrong-way thrusting, and natural hazards in Cascadia. Geology, 29(5): 379-382. DOI:10.1130/0091-7613(2001)029<0379:NCWWWT>2.0.CO;2
Hardy S, Finch E. 2005. Discrete element modelling of detachment folding. Basin Research, 17(4): 507-520. DOI:10.1111/bre.2005.17.issue-4
Hubbert M K, Rubey W W. 1959. Role of fluid pressure in mechanics of overthrust faulting:Ⅰ. Mechanics of fluid-filled porous solids and its application to overthrust faulting. GSA Bulletin, 70(2): 115-166. DOI:10.1130/0016-7606(1959)70[115:ROFPIM]2.0.CO;2
Li W H, Gao R, W ang H Y, et al. 2017. Crustal structure beneath the Liupanshan fault zone and adjacent regions. Chinese Journal of Geophysics (in Chinese), 60(6): 2265-2278. DOI:10.6038/cjg20170619
Liu H Q, McClay K R, Powell D. 1992. Physical models of thrust wedges.//McClay K R ed. Thrust Tectonics. Dordrecht: Springer, 71-81.
Lujan M, Storti F, Balanyá J C, et al. 2003. Role of décollement material with different rheological properties in the structure of the Aljibe thrust imbricate (Flysch Trough, Gibraltar Arc):an analogue modelling approach. Journal of Structural Geology, 25(6): 867-881. DOI:10.1016/S0191-8141(02)00087-1
Malavieille J. 2010. Impact of erosion, sedimentation, and structural heritage on the structure and kinematics of orogenic wedges:Analog models and case studies. GSA Today, 20(1): 4-10.
Marques F O, Cobbold P R. 2002. Topography as a major factor in the development of arcuate thrust belts:insights from sandbox experiments. Tectonophysics, 348(4): 247-268. DOI:10.1016/S0040-1951(02)00077-X
Mary B C L, Maillot B, Leroy Y M. 2013. Predicting orogenic wedge styles as a function of analogue erosion law and material softening. Geochemistry, Geophysics, Geosystems, 14(10): 4523-4543. DOI:10.1002/ggge.20262
Mugnier J L, Baby P, Colletta B, et al. 1997. Thrust geometry controlled by erosion and sedimentation:a view from analogue models. Geology, 25(5): 427-430. DOI:10.1130/0091-7613(1997)025<0427:TGCBEA>2.3.CO;2
Panian J, Wiltschko D. 2004. Ramp initiation in a thrust wedge. Nature, 427(6975): 624-627. DOI:10.1038/nature02334
Panian J, Wiltschko D. 2007. Ramp initiation and spacing in a homogeneous thrust wedge. Journal of Geophysical Research:Solid Earth, 112(B5): B05417. DOI:10.1029/2004JB003596
Price R A, Mountjoy E W. 1971. Geologic structure of the canadian rocky mountains between bow and athabasca rivers-a progress report.//Wheeler J O ed. Structure of the Southern Canadian Cordillera. Toronto: Geological Association of Canada, Special Paper, 7-26.
Ruh J B, Kaus B J P, Burg J P. 2012. Numerical investigation of deformation mechanics in fold-and-thrust belts:Influence of rheology of single and multiple décollements. Tectonics, 31(3): TC3005. DOI:10.1029/2011TC003047
Ruh J B, Gerya T, Burg J P. 2013. High-resolution 3D numerical modeling of thrust wedges:Influence of décollement strength on transfer zones. Geochemistry, Geophysics, Geosystems, 14(4): 1131-1155. DOI:10.1002/ggge.20085
Shen X Z, Zhou Y Z, Zhang Y S, et al. 2013. Geodynamic significance of the crust structure beneath the northeastern margin of Tibet. Progress in Geophysics (in Chinese), 28(5): 2273-2282. DOI:10.6038/pg20130509
Smit J H W, Brun J P, Sokoutis D. 2003. Deformation of brittle-ductile thrust wedges in experiments and nature. Journal of Geophysical Research, 108(B10): 2480. DOI:10.1029/2002JB002190
Storti F, McClay K. 1995. Influence of syntectonic sedimentation on thrust wedges in analogue models. Geology, 23(11): 999-1002. DOI:10.1130/0091-7613(1995)023<0999:IOSSOT>2.3.CO;2
Storti F, Salvini F, McClay K. 2000. Synchronous and velocity-partitioned thrusting and thrust polarity reversal in experimentally produced, doubly-vergent thrust wedges:Implications for natural orogens. Tectonics, 19(2): 378-396. DOI:10.1029/1998TC001079
Wang H Y, Gao R, Yin A, et al. 2012. Deep structure geometry features of Haiyuan Fault and deformation of the crust revealed by deep seismic reflection profiling. Chinese Journal of Geophysics (in Chinese), 55(12): 3902-3909. DOI:10.6038/j.issn.0001-5733.2012.12.003
Wang P, Cui J Z. 2006. Mechanics of Metal Plastic Forming (in Chinese). Beijing: Metallurgical Industry Press: 60-62.
Wang W E, Kirby K, Zhang P Z, et al. 2013. Tertiary basin evolution along the northeastern margin of the Tibetan Plateau:Evidence for basin formation during Oligocene transtension. GSA Bulletin, 125(3-4): 377-400. DOI:10.1130/B30611.1
Zhang Z, Zhang H, Wang L S, et al. 2018. Concurrent deformation in the Longmen Shan and the Sichuan Basin:A critical wedge captured by limit analysis. Tectonics, 37(1): 283-304. DOI:10.1002/2017TC004791
Zhou J X, Xu F Y, Wei C G, et al. 2007. Shortening of analogue models with contractive substrata:Insights into the origin of purely landward-vergent thrusting wedge along the Cascadia subduction zone and the deformation evolution of Himalayan-Tibetan orogen. Earth and Planetary Science Letters, 260(1-2): 313-327. DOI:10.1016/j.jpgl.2007.05.048
蔡永恩, 罗纲, 刘今朝, 等. 2004. 不连续岩体-拱坝系统的动力过程数值模拟方法. 岩土力学, 25(S2): 361-365.
陈九辉, 刘启元, 李顺成, 等. 2005. 青藏高原东北缘一鄂尔多斯地块地壳上地幔S波速度结构. 地球物理学报, 48(2): 333-342. DOI:10.3321/j.issn:0001-5733.2005.02.015
郭晓玉, 高锐, 高建荣, 等. 2017. 综合数据分析青藏高原东北缘六盘山地区构造形变及其构造成因独特性探讨. 地球物理学报, 60(6): 2058-2067. DOI:10.6038/cjg20170603
李文辉, 高锐, 王海燕, 等. 2017. 六盘山断裂带及其邻区地壳结构. 地球物理学报, 60(6): 2265-2278. DOI:10.6038/cjg20170619
沈旭章, 周元泽, 张元生, 等. 2013. 青藏高原东北缘地壳结构变化的地球动力学意义. 地球物理学进展, 28(5): 2273-2282. DOI:10.6038/pg20130509
王海燕, 高锐, 尹安, 等. 2012. 深地震反射剖面揭示的海原断裂带深部几何形态与地壳形变. 地球物理学报, 55(12): 3902-3909. DOI:10.6038/j.issn.0001-5733.2012.12.003
王平, 崔建忠. 2006. 金属塑性成形力学. 北京: 冶金工业出版社: 60-62.